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Preface 


The  original  intent  of  this  effort  was  to  produce  a  fractional 
order  sensor  for  use  in  monitoring  structural  vibrations.  As  the  work 
progressed  it  was  expanded  to  include  analog  simulation  and  proof-oi- 
concept  for  fractional-order  feedback.  This  goal  was  achieved. 
Incorporating  fractional-order  circuits  into  active  feedback  control 
systems  will  provide  many  advantages  for  controls  engineers. 

Many  thanks  to  my  advisor,  Lt  Col  Ron  Bagley.  I  don’t  know  if  I 
can  ever  thank  you  for  your  patience,  understanding  and  commitment. 

Dr.  Torvik,  thank  you  for  insulating  me  from  the  repercussions  that 
could  have  been  there  for  not  finishing  on  time.  Thank  you,  Maj  Kolesar 
for  your  comments  on  the  A/D  recorder.  It  was  the  key  to  our  success. 

To  members  of  my  family  (the  Fischer  clan)  and  Emmanuel  Lutheran 
Church,  thanks  for  your  prayers  and  support.  Special  thanks  to  Wil 
Schonscheck  for  weekly  putting  life  into  proper  perspective. 

To  my  wife  and  kids:  Nita,  your  patience,  encouragement  and  love 
were  there  when  I  needed  them  most.  Nathan  and  Katie,  teenagers 
deserve  more  than  a  dad  whose  nose  is  buried  in  a  book  except  when 
telling  you  to  clean  your  room  and  wash  dishes.  Josh  and  Gabe  -  Daddy 
will  finally  be  able  to  play  soccer  with  you  again. 

Finally,  the  real  credit  belongs  to  Jesus  Christ  who’s  given  me  all 
I  have  and  sees  me  through  the  trials  -  even  though  at  times  the 
battles  seem  overwhelming.  Through  Him,  the  final  victory’s  ours. 
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Abstract 


SUMMARY:  Feedback  of  the  1/2  and  3/2  derivatives  as  veil  as  X  and 

x  is  demonstrated  for  a  second-order  system  defined  by  the  differential 


»quation: 


3/2  .  1/2 

mx  +  c  D  <x>  +c  x  +  c  D  <x>  +  kx  =  u<t> 
It  2  3  t 


Three  methods  of  producing  the  fractional  derivative  or  integral 
of  an  input  signal  are  investigated.  The  method  selected  employs  a 
circuit  developed  at  Trent  University,  Ontario,  Canada  for  use  in 
electrochemistry  research The  circuit  performs  the  mathematical 


i>  v>  / 

operation  d  {JJji*- -"for  -1  <  v  <  1;  negative  values  of  v  represent 
integrationS- The  results  presented  show  the  circuit  accurately 
differentiates  a  sinusoidal  input  for  a  frequency  range  spanning  0.01  Hz 
to  lO.OJlz^) 

^  The  second-order  differential  equation  above  is  simulated  on 
an  analog  computer.  An  optimal  u(t)  is  then  used  for  feedback 
modification  of  the  original  open-loop  system.  Improved  system 
performance  resulted. 


A  Laplace  transform  method  and  a  Mittag-Lef f ler  expansion  provide 
analytical  predictions  of  the  system’s  response.  The  output  of  the  two 
methods  is  identical.  Comparison  of  the  theoretical  predictions  with 
the  experimental  data  shows  excellent  agreement  with  respect  to  the 
initial  transient  behavior  and  asymptotic  behavior  of  the  steady-state 
response  for  both  the  open-  and  closed-loop  systems.  lfrr< 


i  -  ’ 


FRACTIONAL-ORDER  FEEDBACK  IN  LINEAR  SYSTEMS 


I.  Introduction 


d^E  f (x)  ]  (1) 


Members  of  the  scientific  and  engineering  community  have  little 
difficulty  in  recognizing  Eq  (1)  as  the  ^th  derivative  of  a  function  of 
the  variable  x  with  respect  to  time.  But  how  many  will  pose  the 
question  L’Hospital  asked  Leibniz  in  1695,  ’What  if  v  be  one-half  ?'. 
Leibniz  replied  that  it  will  'lead  to  a  paradox',  but  added,  'someday  it 
would  lead  to  useful  consequences'  (1:115).  (For  an  historical 
development  of  fractional  calculus  the  reader  is  directed  to 
reference  1.) 

Fractional  derivatives  have  indeed  proven  useful  in  the  analysis  of 
a  wide  variety  of  physical  systems  (2:126).  In  the  context  of  this 
research,  fractional  derivatives  are  capable  of  describing  the  behavior 
of  viscoelastical ly-damped  structures  (3:348,4:209).  The  spacecraft  and 


(  1  ) 


in  their  designs.  This  expanding  interest  has  motivated  an 
investigation  of  active  control  requirements  for  viscoelastical ly-damped 
systems.  The  stability  of  feedback  in  such  systems  has  been  addressed 
(3:351).  However,  implementation  of  such  a  control  system  requires 
development  of  a  feedback  method  (6)  and  a  sensor  or  instrument  capable 
of  producing  the  fractional  derivative  or  integral  of  electrical 
signals . 

Problem  Statement 

Current  methods  of  providing  the  fractional  derivative  of  an  input 
signal  require  either  digital  processing  or  resistor-capacitor  circuit 
ladders  containing  many  components.  In  addition,  the  resistor-capacitor 
circuit  ladder  is  designed  for  a  fractional  order  of  one-half  (7:39). 

The  goal  for  this  research  was  to  produce  a  general  fractional 
differentiator/integrator  which  serves  as  the  key  component  in  an  active 
control  system  for  viscoelastical ly-damped  structures. 

Approach 

Three  parallel  areas  of  investigation  were  examined  in  an  attempt 
to  identify  a  device: 

1)  Excess  charge  carriers  in  doped  semiconductors  undergo  a 
diffusion  process  in  passing  through  the  semiconductor  lattice.  Bagley 
and  Torvik  showed  fractional  derivatives  of  order  one-half  are  common  to 
diffusion  processes  (8:2).  The  application  of  a  focused  monochromatic 
light  beam  generates  a  localized  excess  concentration  of  charge  carriers 


(  2  ) 


in  the  semiconductor  material  and  results  in  voltage  changes  as  the 
excess  carriers  diffuse  away  from  the  source  (9:54).  Measurement  of 
this  voltage  might  provide  a  possible  source  of  a  fractional  derivative. 

2)  Oldham  identified  a  fractional  integrating  circuit  of  order 
one-half  in  his  textbook  on  fractional  calculus  (10.149).  Thi3  circuit 
has  several  characteristics  which  preclude  it  as  an  optimum  choice. 
However,  in  developing  the  theory  for  this  circuit,  reference  was  made 
to  Wall’s  work  in  continued  fractions  which  can  be  used  to  represent 
Laplace  transforms  of  functions  (11:355).  This  method  might  permit  the 
development  of  ladder  circuits  of  fractional  order,  since  sa  is  the 
Laplace  transform  definition  of  a  fractional  derivative. 

3)  Finally  a  continued  investigation  of  the  literature  to  identify 
other  scientific  fields  that  might  have  fabricated  or  used  such  a  device 
for  research  will  be  undertaken.  Electrochemistry  and  geophysics  are 
two  fields  in  which  fractional  derivatives  have  found  use 
(2:126,10:154) . 

Scope 

The  desirable  attributes  of  the  device  include:  i)  Operation  over  a 
wide  range  of  frequencies  (especially  .01  to  200  Hz  -  the  frequencies  of 
interest  to  many  structural  engineers),  ii)  small  size,  relatively  few 
components,  passive  in  nature  and  inexpensive  to  manufacture,  and 
iii)  generalization  to  any  fractional  order. 


Evaluation  of  device  performance  will  be  accomplished  in  three 
phases : 

Phase  1.  A  function  generator  is  used  to  apply  a  sinusoidal  input 
signal  of  varying  frequency  to  the  device.  The  observed  output  should 
match  the  expected  magnitude  and  phase  behavior  of  a  l/nth-order 
fractional  dif ferent iator/ integrator . 

Phase  2.  An  analog  computer  simulation  of  a  second-order  system 
incorporates  the  device  to  supply  the  fractional  derivative  terms  of  the 
system.  Solution  techniques  for  this  type  of  system  exist  (12,13:141), 
and  experimental  agreement  within  10%  of  the  predicted  response  is 
established  as  the  success  criteria  for  this  phase. 

Phase  3.  This  thesis  parallels  the  development  of  optimal 
fractional-order  feedback  theory  by  Capt  Rich  Walker  (6) .  The  parameters 
developed  for  feedback  by  application  of  his  modified  linear  quadratic 
regulator  will  be  applied  to  the  system  referred  to  in  Phase  2.  The 
ability  of  the  feedback  method  to  produce  improved,  stable  control 
system  performance  will  be  the  criterion  for  success  in  this  phase. 
Assumptions 

This  research  takes  on  only  proof-of-concept  development. 

Circuits,  if  built,  will  consist  of  discrete  components  fabricated  on 
breadboards.  No  attempt  to  use  integrated  circuit  resistors  and 
capacitors  will  be  made.  Perfection  of  circuits  and  device  operation 
will  be  left  as  a  follow-on  task.  Feedback  will  be  accomplished  through 
manual  adjustment  of  attenuators  on  the  analog  computer  and  not  through 
automatic  means.  Most  systems  exhibiting  behavior  that  can  be 


(  4  ) 


modeled  with  fractional  derivatives  involve  the  half-order  derivatives. 
Initial  investigation  will  consist  of  half-order  devices  that  can  be 
generalized  to  l/nth-order  devices.  While  many  materials  exhibit 
stress-strain  and  damping  characteristics  that  can  be  modeled  by 
fractional  derivatives,  no  effort  will  be  made  to  tie  the  system 
simulated  on  the  analog  computer  to  a  physical  system,  or  specific 
material . 


II •  Background 


Evolution  of  Fractional  Calculus 

Beginning  in  1695,  when  L’Hospital  first  posed  the  fractional 
derivative  question,  many  noted  scientists  and  mathematicians  have 
focused  their  attention  toward  developing  the  foundation  of  fractional 
calculus.  In  1823  Neils  Henrik  Abel  used  fractional  calculus  to 
formulate  a  solution  to  the  tautochrone  problem.  This  was  the  first 
recorded  application  of  fractional  calculus  (1:121).  Beginning  in  1832 
Joseph  Liouville  published  several  papers  dealing  with  fractional 
derivatives  (1:116).  G.F.  Bernhard  Riemann  developed  a  theory  for 
fractional  integration  which  was  published  posthumously  in  1876  (1:116). 
But  neither  Riemann  nor  Liouville  were  able  to  provide  definitions  for 
Eq  (1)  which  could  be  generalized  for  any  V,  positive  or  negative,  and 
for  a  sufficiently  wide  class  of  functions.  Henri  Laurent  published 
what  many  consider  the  foundation  paper  in  fractional  calculus  in  1884. 
In  this  paper,  Laurent  produced  a  definition  of  fractional  operations 
which  also  applied  to  integer  values  of  v  and  followed  precisely  the 
rules  of  integer  calculus  for  differentiation  and  integration  (1:118). 
Finally  in  1936,  with  the  theory  of  fractional  operations  well  defined, 
Harold  T.  Davis  developed  a  notation  for  fractional  operations  which 
will  be  used  throughout  the  remainder  of  this  paper  (1:117): 

-v 

Let  V  be  a  positive  real  number,  cDx  f<x>  vill  denote 
integration  of  order  V  of  function  f  along  the  x  axis  from  c 
to  x.  Similarly,  the  operator  cDx  f <x>  denotee  differentia¬ 
tion  of  order  V  with  terminal  limit*  c  and  x.<i:ii7> 


(6) 


The  mathematical  interpretation  of  the  integration  operator  in  the 
time  domain: 

i  t 

D  ~V{i  (t)  3  =  -  f  (t-r)  v~Xi  (r)dr  (2) 

c  T(v)  c 

One  would  assume  that  the  definition  of  fractional  differentiation 
would  involve  replacing  v  by  -v  in  Eq  (2).  However,  this  produces  a 
divergent  integral.  To  obtain  a  fractional  derivative,  first  integrate 
to  the  correct  fractional  order  between  0  and  1  and  then  use 
conventional  differentiation  to  obtain  the  desired  result: 

■  -7 V  131 

d  t  *■  J 

Where  m  is  the  least  integer  greater  than  v,  v  -  m-p,  and  0  <  p  ^  1 
(1:116) . 

Except  for  a  small  class  of  functions,  Eqs  (2)  and  (3)  are  not 
computationally  convenient  to  use.  However,  if  both  expressions  are 
treated  as  convolution  integrals,  the  fractional  operators  have  simple 
definitions  in  the  Laplace  domain: 

£{  cDt_l,f(t)  >  =  a'VX(  f(t)  >  (4) 

£{  cDtVf(t)  )  =  sm-p*  f(t)  )  (5) 

It  is  through  the  use  of  these  Laplace  transform  operators  that 
fractional  calculus  proves  to  be  most  useful  (2.126). 
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Abel  used  fractional  calculus  to  determine  the  shape  of  a  wire  in  a 
vertical  plane  which  would  cause  a  bead  placed  on  the  wire  to  take  the 
same  amount  of  time  to  reach  the  lower  end  no  matter  where  the  bead  was 
originally  placed  (1:121).  In  1921  P.G.  Nutting  noted  stress  relaxation 
properties  in  viscoelastic  materials  could  be  modeled  by  fractional 
powers  of  time  rather  than  the  traditionally  expected  decaying 
exponential  (2:126).  This  development  led  A.  Gemant  in  1938  to  propose 
time  differentiation  of  fractional  order  for  modeling  stiffness  and 
damping  properties  of  viscoelastic  materials  (2:126).  In  1966, 

M.  Caputo  suggested  fractional  derivatives  might  be  used  to  model  the 
behavior  of  geological  strata  (2:126).  In  1970,  V.  0.  Shestopol 
employed  fractional  calculus  to  describe  the  mechanism  of  deformation  in 
tungsten  and  platinum  at  high  temperatures  (2:126).  Other  areas  of 
physical  application  include  creep  and  stress  relaxation,  creep 
buckling,  and  techniques  for  fitting  experimental  data.  (2:126  -  This 
reference  is  an  excellent  source  for  a  historical  perspective  of  the 
application  of  fractional  calculus  to  damping.) 

Beginning  in  1979,  R.  Bagley  and  P.  Torvik  published  a  series  of 
articles  developing  a  method  for  modeling  the  behavior  of  viscoelastic 
materials  by  using  fractional  derivatives  in  a  finite  element 
formulation.  The  stress-strain  curve  for  viscoelastic  materials  is 
dependent,  not  only  on  the  applied  load,  but  also  on  the  frequency  at 
which  the  load  is  applied.  Previous  methods  of  modeling  this  phenomena 
were  computationally  difficult  to  use  for  anything  but  steady-state 
conditions  or  predicted  non-causal  behavior  --  both  of  which  made 


t  transient  analysis  of  viscoelastic  systems  impractical  (14.741).  Bagley 

showed  that  a  fractional  derivative  model,  incorporating  three  to  five 
parameters,  could  be  used  to  accurately  describe  the  behavior  of  many 
i  viscoelastic  materials  over  a  large  range  of  frequencies  (2:128). 

Approximately  130  materials  have  been  characterized  in  this  manner 
(2:130).  Subsequent  research  verified  the  damping  characteristics  of 
i  viscoelastic  materials  could  also  be  modeled  in  a  similar  fashion 

(15:83).  In  1985,  Bagley  and  Torvik  showed  finite  element  methods  could 
incorporate  fractional  calculus,  effectively  increasing  the  usefulness 
of  fractional  calculus  to  the  structural  analyst  (14:743). 

When  a  model  effectively  describes  a  physical  system  over  an  entire 
range  of  operating  conditions,  one  might  expect  the  model  to  also 
predict  behavior  for  conditions  outside  of  range  tested.  One  of  the 
areas  of  recent  interest  in  engineering  is  the  damping  requirements  of 
large  space  structures.  The  excellent  damping  properties  of 
viscoelastic  materials  make  them  highly  suitable  candidates  for  use  in 
such  structures.  The  space  environment  cannot  be  simulated  in 
ground-based  testing.  A  model  that  can  predict  performance  outside  of 
testable  conditions  is  necessary.  Bagley  and  Torvik  present  a  strong 
argument  for  accepting  the  fractional  calculus  model  as  such  (2:134). 

Requirements  for  these  structures  will  no  doubt  include  some 
combination  of  active  and  passive  control  for  station-keeping  and  other 
mission  requirements.  If  viscoelastic  materials  are  integral  to  these 
structures,  it  seems  likely  that  the  control  system  should  include  and 
adapt  the  fractional-order  model.  Three  questions  immediately  arise  -- 
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can  fractional  orders  of  displacement  be  sensed;  is  there  a  constitutive 
control  law  for  fractional-order  feedback  in  these  systems;  and,  is  such 
feedback  stable  “> 

Steven  B  Skaar,  et  al ,  addressed  the  stability  of  fractional  order 
feedback  in  simple  physical  systems  and  outlined  conditions  necessary  to 
guarantee  a  stable  system  (3) .  The  paper  discussed  the  characteristics 
of  a  fractional-order  sensor  and  suggested  a  circuit  comprised  of 
lead-lag  transfer  functions  (3:355).  However,  the  authors  stated  more 
study  was  required  to  develop  the  idea.  In  addition,  the  paper 
addressed  the  form  of  the  feedback  control  law.  The  feedback  must  be  a 
function  of  the  vibration  frequency  and  must  be  of  the  same  order  as  the 
equations  of  motion.  In  other  words,  if  the  material  properties  are 
best  described  by  the  one-half  and  three-halves  derivatives,  the 
feedback  should  be  of  the  same  form  (3:356). 

Capt  Rich  Walker  went  further  in  his  1988  AFIT  Master's  Thesis, 
defining  an  optimal  linear  quadratic  feedback  control  law.  His  research 
showed  that  linear  feedback  of  the  fractional-  and  integer-order  terms 
of  a  second-order  system  could  both  improve  performance  and  guarantee 
stability.  One  of  the  assumptions  made  in  his  paper  was  the  existence 
of  a  sensor  or  device  which  would  provide  the  fractional-order  input 
signals  to  the  controller  (6). 

Summation 

Models  exist,  using  fractional  derivatives,  which  accurately 
describe  and  predict  the  behavior  of  viscoelastic  systems.  The  use  of 
viscoelastic  systems  in  large  space  structures  would  provide  a  number  of 
benefits.  A  control  law  exists  for  feedback  control  of  systems 
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containing  these  viscoelastic  materials.  This  law  improves  system 
performance  and  produces  stable  systems.  To  implement  this  control  law, 
a  device  is  required  which  produces  the  fractional-order  derivative  of 
an  input  signal.  This  is  the  justification  of  the  work  to  follow  in 
this  thesis. 

Potential  Sources  of  Development 

A  rigid  flat  plate  on  the  surface  of  a  Newtonian  fluid  produces 
motion  in  the  fluid  which  can  be  described  by  the  diffusion  equation 
(8:2)  : 


<?t 
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v  sb  transverse  velocity  profile  in  the  fluid 
z  =  vertical  distance  from  the  surface  of  the  plate 
a  =  ratio  of  fluid  viscosity  to  fluid  density 
t  =  time 

Simply  stated,  the  time  rate  of  change  of  the  fluid  is  proportional 
to  the  second  derivative  of  velocity  with  respect  to  distance  from  the 
surface  of  the  plate.  Bagley  and  Torvik  show  that  the  solution  of  this 
equation  is  a  fractional  derivative  of  order  one-half  (8:2).  Similar 
solutions  can  be  found  for  other  parabolic  differential  equations.  One 
would  expect  that  any  diffusion  process  would  result  in  fractional-order 
solutions . 

The  goal  is  to  find  a  diffusion  process  where  current,  voltage,  or 


charge  carriers,  excited  by  an  electrical  or  mechanical  signal,  diffuse 


through  some  medium.  The  Haynes-Shockley  experiment,  circa  1949,  is 
used  to  create  excess  charge  carriers  in  semiconductor  materials  by  the 
application  of  a  focused  monochromatic  light  source.  These  excess 
charge  carriers  then  diffuse  through  the  semiconductor  lattice  and  can 
be  measured  at  some  distance  from  the  source  of  the  excitation  (9:56). 
One  might  expect  the  charge  measured  to  be  related  to  the  input  source 
by  some  type  of  diffusion  process.  This  will  be  one  area  of 
investigation . 

Another  type  of  diffusion  process  provides  a  second  possible  device 
candidate.  Electrochemists  analyze  solutions  for  electro-reducible 
species  with  a  two-electrode  cell  configuration.  When  the  potential 
(voltage)  of  one  of  the  electrodes  is  lowered,  electro-reduction  occurs 
at  this  electrode  which  creates  a  time-dependent  faradaic  current.  The 
half-order  integral  (semi integral )  of  this  current  is  proportional  to 
the  bulk  concentration  of  the  solution’s  oxidizer.  The  first  method 
used  to  obtain  the  semi  integral  involved  recording  the  current  versus 
time  response  of  the  system;  the  result  was  then  digitally  processed. 

K.  Oldham  developed  a  resistor/capacitor  ladder  circuit  that 
produced  the  analog  semi  integral  when  hard-wired  into  the  reduction 
apparatus.  This  innovation  significantly  reduced  the  sample  processing 
time  of  the  analysis.  Oldham  states  this  technique  has  now  become  the 
preferred  method  (7:39).  One  difficulty  with  Oldham’s  circuit  is  the 
number  of  cells  (resistor/capacitor  combinations)  required  to  produce  an 
accurate  output.  Another  is  its  limitation  to  half-order  fractional 
operations,  instead  of  general  fractional  orders. 
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The  theory  of  continued  fractions  developed  by  Wall  (11)  was 
significant  in  Oldham's  development.  Use  of  this  theory  might  possibly 
lead  to  the  synthesis  of  a  passive  component  circuit  which  will  produce 
a  generalized  fractional-order  output.  This  will  be  the  alternative 
area  of  investigation. 
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Ill .  Theory 


A  specialized  form  of  the  continuity  equation  for  excess  minority 
carriers  in  semiconductors  is  (9:56): 


#n  n-n 

_ P  _  _  po 

9t  T 
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_ p 

at 
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where  np  =  excess  minority  carrier  density  (free  electrons)  in 
p-type  semiconductor 

n?o  s  excess  minority  carrier  density  at  equilibrium 

M  =  minority  carrier  mobility 
n 

£  s  applied  electric  field 

s  minority  carrier  diffusion  coefficient 
r  =  minority  carrier  lifetime. 

n 

This  equation  arises  from  pulsing  the  semiconductor  surface  with  focused 
beam  of  monochromatic  light  of  the  proper  wavelength  to  generate  excess 
carriers.  The  term  (np-np0)/T^  rePresents  the  generation  rate  per 
unit  volume  due  to  the  excitation  pulses.  The  carrier  lifetime,  t 
represents  the  average  time  a  minority  carrier  exists  before 
recombining  with  a  hole.  This  parameter  is  a  direct  measure  of  the 
number  of  impurities  and  imperfections  in  the  semiconductor  lattice. 

The  greater  the  number  of  imperfections,  the  shorter  the  lifetime. 
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The  solution  to  Eq  (8)  : 


n  (x,t) 


N 


(47TD  t)1/Z 

n 


(9) 


where  N  =  number  of  carriers  generated  per  unit  volume  (9:55). 

1/2 

The  term  N/(4rtD  t)  is  in  the  characteristic  form  of  the  result 

n 

associated  with  a  half-order  derivative  operation.  However,  the 
exponential  term  multiplying  it  gives  a  decaying  exponential  appearance 
to  the  excess  minority  carrier  population.  If  the  exponents  could  be 
set  approximately  equal  to  zero,  the  negative  voltage  from  the  free 
electrons  would  be  inversely  proportional  to  the  half-order  derivative 
of  the  number  of  free  electrons  generated  per  unit  volume.  In  an  analog 
simulation,  the  strength  of  the  input  pulse  would  be  proportional  to  the 
voltage  of  the  input  signal.  The  rate  at  which  the  carriers  decay  at 
the  measurement  point,  x,  would  be  related  to  the  half-order  derivative 
of  the  input.  This  criteria  requires  that  two  conditions  be  satisfied 
since  all  the  terms  in  the  exponent  are  positive.  That  is, 


x  «  0  (10) 

t  =  0.001T 

n 

The  first  condition  requires  that  the  charge  measurement  be  made 
close  to  the  excitation  site.  Inserting  a  probe  close  to  the  excitation 
site  would  affect  the  semiconductor  properties,  invalidating  Eq  (6). 

-9 

Typical  lifetimes  in  silicon  and  germanium  are  on  the  order  of  10 

seconds  (9:849).  Thus,  the  maximum  effective  time  in  which  to  accomplish 

-<$ 

the  measurement  would  be  10  seconds.  The  lowest  operating  frequency 
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of  the  device  would  then  be  in  the  Megahertz  range.  Most  structural 
engineers  are  interested  in  frequencies  three  to  four  orders  of 
magnitude  smaller  than  this.  Some  other  method  is  needed;  this 


diffusion  process  is  not  a  suitable  candidate. 

The  ladder  circuit,  illustrated  in  Figure  1,  which  consists  of 
T-cells  composed  of  resistors  and  capacitors,  closely  approximates  a 
half-order  integrator  if  the  number  of  cells  are  on  the  order  of 
100  (10:154).  Several  related  articles  also  simulate  half-order 
operations  with  simpler  circuitry  (7:41,16:253).  These  other  methods 
were  also  limited  to  half-order  operations. 


figure  1.  T-Cell  Ladder  Circuit- 


One  reference  detailed  the  representation  of  certain  functions 
using  continued  fractions  (11).  One  function  was  of  particular 
interest : 
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1 

1  +  az 
1  +  z 

1  +  <a+l>z 

i  +  2Z  . 

The  expression  on  the  right  hand-side  of  Eq  (11)  is  called  a  continued 

fraction  (11:349).  If  the  following  substitutions  are  made: 

u  =  St  (12) 

z  =  (RCs)  1 
du  =  sdt 


oo  -u . 
r  e  du 


o  , ,  .a 
(1+  zu) 


the  term  on  the  left  becomes: 


s 


(13) 


This  form  looks  like  a  Laplace  transform.  If  the  substitution 
r  =  RC  +  t  is  made,  the  following  transform  results: 


(RC)°s 


RC»  0»  _ 

e  (BC)  Hl-oO 


a 

s 


(14) 


Cl 

From  Section  II,  s  is  recognized  as  the  Laplace  transform  of  the 
fractional  derivative  of  order  a. 
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The  continued  fraction,  with  (RCs)"1  substituted  for  z.  is  the 
impedance  expression  for  the  circuit  illustrated  in  Figure  1.  For  a 
given  circuit,  the  impedance  relation  is  written: 


c  (s) 
i  (s) 


z  (s) 


(15) 


As  an  example,  Figure  2  illustrates  an  abbreviated  version  of  the  ladder 
circuit,  and  the  impedance  relation  becomes: 


z  (s) 


R  + 

l 


(16) 


Figure  2.  Abbreviated  Ladder  Circuit 
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by  inverting  both  aides  of  Eq  (17)  and  substituting  Eq  (15)  for  z(s): 


i(s)Bi  = 

£(s) 


i  + 


R  C  a 

1  O 


1  t 


R  C  a 
O  O 


eRC“(RC)a  Hi-oO 


a 
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But  this  is  true  only  if  the  resistors  and  capacitors  in  the  circuit 
have  values  such  that: 


B  .  =  0=^  -L  B  ;  C  .  -Xiliai.  c  ll9) 

T(l+a)  m!  r(m+a) 

where  n  =  the  number  of  cells  in  the  circuit  and  m  goes  from  0  to  n 
This  finally  produces  an  expression: 

i(s)  R  =  e  (RC)  r(l-a)  s  £(s)  (20) 

This  expression  states  the  voltage  drop  across  the  first  resistor 
is  the  fractional  derivative  of  the  input  voltage  multiplied  by  a 

RGs 

constant.  However,  the  term  e  is  a  delay  term.  If  the  substitution 
is  made  to  transform  s  into  the  frequency  domain,  then  for 

s  =  joo  (21) 

the  delay  term  becomes  a  function  of  frequency,  implying  the  phase 
shift  of  the  circuit  is  not  constant,  but  also  is  a  function  of 
frequency.  This  circuit  was  fabricated  and  tested.  The  Bode  gain  plot 
had  the  proper  gain  slope,  but  the  circuit  had  a  phase  shift  of  45  ±  5 
degrees  over  less  than  a  decade  of  frequency. 
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Further  investigation  of  the  literature  provided  an  article 
describing  a  circuit  which  produced  fractional  derivative  (or 
integral)  output  (17).  The  authors  claimed  the  circuit  operated 
successfully  in  the  millisecond  to  several  hundreds  of  seconds  range. 
Since  frequency  is  the  inverse  of  time,  this  performance 
characteristic  translated  into  frequencies  spanning  a  range  of  0.01  to 
1000  Hz.  However,  the  circuit  had  only  been  used  with  constant  and  ramp 
input  signals.  One  of  the  assumptions  used  in  designing  the  circuit  was 
that  the  input  current  was  constant.  Because  the  circuit  is  a  linear 
device,  it  should  be  useful  in  processing  sinusoidal  signals.  Again, 
the  circuit  contains  chains  of  resistors  and  capacitors,  but  this  time 
the  cells  are  connected  as  illustrated  in  Figure  3.  The  following 
narrative  provides  the  rationale  for  why  the  circuit  works: 


I  f  an  input  current  is  polled  through  th®  network  a  potential 
E  <t>  is  generated  across  th®  jth  pair  of  componsnls.  The 
magnitude  of  th®  current  i®  r®lat®d  to  lh®  component  valu®s 

of  R  and  C.  by: 
j  J 


E .  <  t> 

J  d  ,  , 

ut>  =  — —  +  c  — -  E  E  <t>] 

R.  j  dt  J 

J 


(22) 


which  can  b®  inv®rt®d  <by  Laplac®  transformation): 

t  (23) 

E  <t>  =  /  ia-T>  ®xp  r  - — -  1  dT 

J  C.  «■  R  C  J 

J  c  j  J 


T  boing  an  integration  variable.  Sine®  i<t>  is  common  to  all 
pairs  it  follows  th®  potential  across  th®  entire  circuit  is: 


E<t)  = 


t  +N 

/  i<t-T>  Z  — —  exp  r  -  1  dT 

j  =  -n  C  t-RCJ 

C  j  j  j 


(24) 
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V  in 


Figure  3.  Oldham-Zoski  Circuit  Schematic 

Adjacent  resistors  are  now  stipulated  to  differ  in  vatu® 
by  a  constant  factor,  as  do  the  capacitor  valuss,  although 
the  progression  ratios  are  not  necessarily  equal.  Thus: 

(25) 

R  =  g  JR  and  C  .  =  O  JC 

j  O  j  O 

where  a  and  g  are  geometric  ratios,  both  >1.  A  parameter  V 
is  defined  relating  a  to  g: 

In  a 

(26) 

1  In  g 

and  Eq  <24>  can  be  recast 

E<t)  =  *°  /Ul"T>  +Z  T°JgJ  exp  [  -TaV]  dT 

1-u  V  j  =-n  RC  L  R  C  j 

c  CTJ  C  C  OO 

O 

(27) 

In  the  limit  as  n  and  N  approach  infinity  and  O  and  g 
approach  unity: 
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and  by  using  th*  Rismann-Uiouvt  l  Is  dsfinUion  of  a  fractional 
integral  as  in  S«cti on  II  (17:27)  : 


E<t>  = 


V> 

ft  cac<U7T>  R  v-1  , 

_ O  d  [  i  <  t>  ] 

r-.  ,  i.-V  ,  y-1 

r <1- V>  ln<a  g>  C  dt 


(29) 


Therefore,  the  input  voltage  is  proportional  to  the  integral  of  the 
current  through  the  circuit.  From  the  principle  of  continuity  the  input 
current  equals  the  output  current.  If  a  resistor  is  connected  in  series 


with  this  circuit,  the  voltage  across  the  resistor,  ER(t)  i3  related  to 
the  input  voltage  E(t) : 


ER(t)  =  i(t)  R  = 


ln<Og>  C 


Tl  c»c<  VTl>  R 


dV  E  <  «•>] 


The  voltage  across  this  resistor  is  the  i/vth  derivative  of  the  input 
voltage.  Additional  details  of  this  circuit  -  henceforth  referred  to 


as  the  Oldham-Zoski  circuit  (OZ)  -  are  contained  in  Appendix  E 
(including  fabrication  and  interface  details). 


The  theory  for  simulating  a  second-order  equation,  including 
half-order  terms  supplied  by  the  OZ  circuit,  is  developed  in  Section  IV 
under  the  subheading  Analog  Simulation.  Basically,  given  a  set  of 
initial  conditions  and  a  set  of  coefficients  for  a  second-order  system, 
the  response,  as  a  function  of  time,  can  be  obtained  on  the  analog 
computer . 
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Finally,  there  is  a  need  to  evaluate  experimental  results  and  be 
confident  in  their  accuracy.  Two  techniques  will  be  used  to  predict  the 
theoretical  results  of  the  open  and  closed  loop  simulations.  First,  a 
Laplace  transform  method  which  will  be  used  to  evaluate  the  residues  and 
a  contour  integral  in  the  s  plane.  Second,  a  Mittag-Lef f ler  expansion 
which  will  be  used  to  expand  the  homogeneous  solution  of  the  differential 
equation  as  an  infinite  series  in  the  time  domain.  Details  of  the 
Laplace  transform  method  are  included  in  Appendix  A.  Appendix  C  details 
the  Mittag-Lef f ler  expansion. 
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IV.  Experimental  Procedure 


Overview 

This  section  details  eight  distinct  tasks: 

1)  Circuit  design  and  fabrication.  The  details  of  realizing  the 
half-order  derivative  circuit  are  presented. 

2)  Circuit  Performance  Evaluation.  A  signal  generator  was  used  to 
produce  a  sinusoidal  signal  of  varying  frequency  as  the  input  to  the 
half -order  circuit.  The  output  of  the  circuit  was  compared  to  the  input 
signal  to  ensure  the  proper  phase  and  magnitude  relationships  existed. 

3)  Circuit  Gain  Adjustment  for  Analog  Simulation.  To  obtain 

a  zero  amplitude  offset  of  the  derivative  signal,  several  operational 
amplifiers  and  a  potentiometer  were  combined  with  the  half-order 
circuit . 

4)  Half-Order  Analog  Performance  Evaluation.  A  harmonic 
oscillator  circuit  was  programmed  on  an  analog  computer.  Its  output  was 
introduced  into  the  circuit  and  op  amp  half-order  combinations  to  ensure 
proper  phase  and  magnitude  performance. 

5)  Fractional-Order  Equation  Synthesis.  Using  the  physical  variable 
method,  a  second-order  differential  equation,  including  half-  and 
three-halves  order  terms,  was  programmed  on  an  analog  computer. 

6)  Open-  and  Closed-Loop  Simulation.  The  performance  results  and 
comparison  with  the  analytic  prediction  are  presented.  As  with  all 
experiments  there  were  some  discrepancies .  The  source  and  resolution  of 
these  discrepancies  is  discussed. 
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7)  Total-Cycle  Simulation.  With  the  system  totally  at  rest,  a  step 


# 


input  was  introduced  into  the  system  to  establish  the  initial  conditions 
tested.  This  task  identified  adjustments  to  correctly  model  the 
initial-value  problem. 

8)  Final  Configuration  and  Performance  Evaluation.  The  lessons 
learned  from  the  total-cycle  simulation  were  applied  to  an  initial-value 
configuration.  The  final  results  are  discussed. 

This  section  provides  a  general  functional  description  of  the 
procedures  used  in  each  task  without  going  into  a  detailed  step-by-step 
description  of  the  equipment  configuration  and  operation.  Learning  to 
operate  the  equipment  correctly  required  a  significant  effort  during 
this  research.  Therefore,  where  appropriate,  references  are  made  to 
Appendix  F  which  contains  a  detailed  description  of  the  equipment 
configuration  and  operation. 

Circuit  Design  and  Fabrication 

For  this  application  a  half-order  differentiating  Oldham-Zoski 
circuit  will  be  designed  and  fabricated.  Appendix  E  details  the  method 
of  determining  the  component  values  listed  in  Table  1.  This  circuit  is 
shown  schematically  in  Figure  3.  Resistors  and  capacitors  within  2X  of 
the  design  values  listed  in  the  table  were  used.  This  is  consistent 
with  the  Oldham-Zoski  article  (17:35).  A  photograph  of  the  completed 
circuits  is  shown  in  Figure  4.  The  as-fabricated  component  values  are 
also  listed  in  Table  1.  These  are  not  of f -the-shel i  values,  but  were 
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obtained  by  serial  and  parallel  combination  of  discrete  components. 

This  was  necessary  as  the  design  values  were  not  commercially  available 
Trial  and  error  combinations  of  components  were  evaluated  until  a 
combination  within  21  of  the  desired  value  was  obtained. 


Figure  4.  01 dham-Zoski  Circuits  Constructed 

on  Breadboard 
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After  completion  of  component  selections,  the  combinations  were 
placed  on  an  E&L  Instruments  Elite  1  circuit  board  and  rechecked  to 
ensure  the  values  matched  the  measured  values.  This  step  is  highly 


recommended  because  several  combinations  were  determined  to  be  outside 
the  tolerance  limits  when  installed.  Two  circuits  were  built,  Circuit  2 
was  designed  and  fabricated  with  twice  the  capacitance  and  half  the 
resistance  of  Circuit  1.  Circuit  1  was  used  to  provide  the  half 
derivative,  and  Circuit  2  provided  the  three-halves  derivative  in  the 
analog  simulations. 

Circuit  Performance  Evaluation 

The  simplest  criteria  to  apply  in  identifying  a  fractional 
derivative  of  order  one-half  is  its  Bode  plot.  A  gain  of  lOdB  per 
decade  and  a  constant  phase  shift  of  45  degrees  characterizes  a 
half-order  derivative  (18:226).  Experimentally  determining  a  circuit's 
Bode  plot  requires  a  sine  wave  generator,  a  recording  or  display  device, 
and  a  counter/timer .  The  method  employed  connects  the  output  of  the 
Tektronix  function  generator  (model  FG506)  to  the  input  of  the  circuit 
under  test.  The  output  of  the  circuit  under  test  is  connected  to  either 
a  Hewlett-Packard  Analog/Digital  recorder  (model  7090A)  or  a  Tektronix 
oscilloscope  (model  SC504) ;  a  Tektronix  counter/timer  (model  DC509) 
is  connected  to  the  output  of  the  function  generator.  The  output  voltage 
is  then  compared  to  the  input  voltage,  and  the  gain  is  calculated. 

The  gain  at  a  given  frequency  (to)  ia  given  by: 


Gain , 


ICO 


constant 


__  .  foul  put  voltage"^ 

20  los,o|— - — -J 

'■input  voltage 


(31) 
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The  phase  shift  is  calculated  by  determining  the  time  interval  (tps) 
between  the  t-axis  crossings  of  the  input  and  output  voltages  and 
comparing  it  to  the  total  period  (T)  of  the  sine  wave  as  illustrated  in 
Figure  5. 


Phase  Shift  (deg) 


(32) 


This  procedure  is  repeated  for  a  number  of  frequencies  in  the 
performance  range  of  the  circuit.  The  data  was  then  plotted,  as 
illustrated  in  Figures  6-9  ,  and  the  slope  of  the  gain  plot  and  the 
magnitude  of  the  phase  plot  are  compared  to  the  lOdB/decade  gain  slope 
and  45  degree  phase  shift  criteria.  A  detailed  equipment  list  and  test 
procedure  is  contained  in  Appendix  F. 

Experimental  values  of  10  ±  1  dB/decade  for  the  gain  and  45  ±  5 
degrees  for  the  phase  shift  were  established  as  the  initial  pass/fail 
criteria  for  the  circuit.  These  values  were  a  first  attempt  at 
establishing  such  a  criteria.  Since  both  circuits  tested  satisfied  this 
criteria  and  performed  adequately  on  the  analog  computer,  these 
parameters  appear  to  be  a  valid  performance  criteria. 
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Figure  6.  Cir 


PHASE  ANGLE  (DEG) 


Circuit  Gain  and  Adjustment  for  Analog  Simulation 

The  change  in  gain  and  phase  of  the  differentiating  circuit’s 
output  arises  from  the  mathematical  definition  of  the  steady-state 
derivative  of  a  sine  function: 

av  v  (33) 

—  (A  sin  cot)  =  Aco  (sin(tot  -  vn/z)  ] 
dt 

When  the  input  amplitude  is  unity,  the  magnification  factor  of  the 
output  amplitude  with  regards  to  the  input  amplitude  is  co  -  which  is  a 
function  of  the  frequency.  If  the  Oldham-Zoski  circuit  is  to  function 
as  a  true  differentiator,  this  amplitude  magnification  factor  at  u>  =  1 
rad/sec  must  be  \V  -  1.  Figure  6  and  Figure  8  show  clearly  that  this 
is  not  the  case  for  either  of  these  circuits.  Therefore,  both  circuits 
must  be  combined  with  operational  amplifiers  to  ensure  that  the  gain 
plot  has  a  value  of  zero  dB  at  w  =  i  radian/sec. 

The  Oldham-Zoski  article  defines  relationships  for  the 
amplification  required  for  a  given  circuit  design  (17:30).  For  Circuit 
1  a  gain  of  26.53  (pure  number,  not  dB)  was  predicted;  the  gain  for 
Circuit  2  was  13.54.  To  realize  this  on  the  analog  computer,  two 
amplifiers  with  gains  of  10  and  a  potentiometer  in  series  with  each 
half-order  circuit  were  required  as  shown  in  Figure  10.  More 
information  on  the  gain  calculations  for  each  circuit  is  contained  in 


Appendix  E. 


Figure  10.  Operational  Amplifier  and  Oldham-Zoski 
Circuit  Combination  for  Q^(x) 


Half  ^  Order  Circuit  Analog  Perf ormance  Evaluation 

Theory  does  not  always  match  real  life,  and  verification  of  the 
gains  was  required  prior  to  simulating  a  second-order  system.  A  sine 
wave  with  w  =  1  radian/sec  was  used  to  perform  the  verification.  A 
harmonic  oscillator  was  programmed  with  a  frequency  of  1  radian/sec 
(see  Figure  11  and  reference  Appendix  F  for  details  regarding  the 
programming  of  a  harmonic  oscillator  on  the  analog  computer)  .  The 
oscillator's  output  served  as  the  input  to  the  half-order  circuit  being 
evaluated.  The  input  amplitude  and  the  output  amplitude  were  plotted 
using  the  Hewlett-Packard  A/D  recorder  (  model  7090A  -  see  Figure  12) .  A 
comparison  is  made,  and  adjustments,  if  necessary,  are  determined  by  the 
relationship: 

(34) 

actual  output  amplitude  i 

input  amplitude  adjustment  factor 
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Figure  11.  Harmonic  Oscillator  -  Analog  Configuration 


Figure  12.  Amplitude  Validation  Circuit  Using 
Harmonic  Oscillator 


The  initial  potentiometer  setting  for  the  half-order  circuit 

was  modified  by  multiplying  its  value  by  the  adjustment  factor.  The 

process  was  repeated  until  the  adjustment  factor  was  approximately  one. 

An  example  of  this  process  is  depicted  in  Figures  13  and  14. 

Although  the  circuits  are  tuned  to  perform  correctly  at  co  =  1 

radian/sec,  their  behavior  at  other  frequencies  requires  evaluation.  As 

2 

a  check,  the  frequency  of  the  harmonic  oscillator  was  changed  to  co  =  3.5 
radians/sec.  From  Eq  (33)  the  predicted  magnification  factor  for  the 
output  amplitude  at  this  frequency  should  be  (3. 5) 1/4  or  1.37.  Figure 
15  shows  the  experimental  result  of  1.35.  The  combination  of  this 
check,  and  the  original  circuit  validation,  bolstered  confidence  in  the 
half  derivative  circuits  prior  to  attempting  the  analog  simulations. 
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Initial  Dy(x)  .  Utlx)  Using  Theoretical  Predictions 
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Figure  14.  D?(x)  ,  O^x)  After  Gain  Adjustments 
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Figure  15.  Dy(x)  .  Ut\x)  Gain  Check  at  uJ  ^  =  3.5  rad/sec 


It  would  be  desirable  if  both  the  1/2  and  s^2  derivatives  could  be 
obtained  by  using  a  single  half-order  circuit.  Using  the  Oldham-Zoski 
circuit  to  differentiate  the  x  signal  to  provide  Dt  (x) ,  and  then 

3/2 

connecting  the  Dt  (x)  output  to  a  full  integer  integrator  should 

produce  the  D*^2(x)  signal.  All  analog  integrators  are  checked  at  the 

factory  to  ensure  their  output  falls  within  a  range  that  matches  the 

performance  of  the  others.  Using  only  one  half-order  circuit  to 

1/2  3/2 

provide  the  signal  for  both  Dt  (x)  and  Dt  (x)  eliminates  concerns  with 
performance  matching  two  half-order  circuits.  To  test  this  concept  the 
harmonic  oscillator  was  tuned  to  <0  =  1  radian/sec,  its  x  signal  was 
differentiated  to  obtain  (x) ,  then  integrated  to  obtain  Dt  (x)  . 

This  circuit  configuration  is  pictured  in  Figure  16,  and  the  results  are 
plotted  in  Figure  17.  The  peak-to-peak  amplitude  for  each  cycle  is 
correct,  and  the  phase  shift  is  also  correct.  However,  there  is  an 
underlying  ramp  in  the  signal,  indicating  a  bias  voltage  exists  in  the 

3/2 

Dt  (x)  signal.  Several  possible  causes  were  identified,  (see  Section 
V)  and  fixes  were  implemented.  Even  then,  the  ramp  could  not  be 
totally  eliminated. 

1/2 

Incorporating  a  second  half-order  differentiator  to  produce  (x) 
from  the  x  (displacement)  signal  solved  much  of  the  problem.  With  the 
system  configured  as  shown  in  Figure  18,  the  D*/2(x) ,  x,  and  D*^2(x) 
signals  are  shown  in  Figure  19.  The  proper  relationships  now  exist. 
Section  V  includes  a  discussion  of  the  problems  associated  with  using 
a  single  circuit  to  produce  both  fractional  derivatives. 
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Figure  16.  Configuration  for  DJtx)  .  Dftx)  Using  Only  One 
Fractional  Order  Circuit 
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Figure  17.  Results  of  Integrating  D^C X )  to  Obtai 
Dtflx)  Directly 
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Fractional-Order  Equation  Synthesis 


8 

The  EAI  analog  computer  (PACE, model  TR  48)  depicted  in  Figure 
20  is  a  solid-state  second-generation  computer.  Components  available 
for  patching  into  simulations  include  operational  amplifiers  capable  of 
summing,  inverting,  multiplying,  and  integrating;  potentiometers  for 
precision  multiplication;  comparators  for  switching  logic;  function 
switches  for  applying  various  input  (forcing)  functions;  and  a  digital 
voltmeter  for  readout  purposes.  Each  component  has  input  and  output 
jacks  which  accept  patch  cords  with  banana  plug  terminations  (see  Figure 
21) ,  permitting  interconnection  of  components  for  modeling  a  variety  of 
systems.  Appendix  F  contains  detailed  descriptions  of  these  components 
and  their  uses. 


Fiaure  20.  Pace  TR48  Analog  Computer 
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Figure  21.  Patch  Cords  and  Input  Jacks 


A  general  second-order  differential  equation  containing  half-order 
derivatives  can  be  written  as  follows: 


mx 


♦  Cjp*/2(x) 


c  *  + 
2 


ft 


+  kx  *  u(t) 


(35) 


where 


m  =  system  mass 
cz  =  viscous  damping  coefficient 
k  =  system  stiffnes  coefficient 
ci  =  viscoinertial  damping  coefficient 
cg  =  viscoelastic  damping  coefficient 
u(t)s  system  control  or  forcing  function 
The  x,  x,  and  x  terms  are  acceleration,  velocity,  and  position.  The  two 

3/2  1/2 

additional  terms,  Dt  (x)  and  Dt  (x) ,  are  respectively,  the 
three-halves  and  one-half  derivative  of  position  with  respect  to  time. 


* 
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Eq  (35)  can  be  rewritten  such  that  the  acceleration  is  a  function 
of  the  other  terms: 

x  =  m  1  £  f  (t)  -  c±D2/2(x)-  c2h  -  cj)^x)  j  (36) 

The  three-halves  derivative  can  be  obtained  by  differentiating  x 
or  integrating  x  by  order  one-half,  and  is  mathematically  defined  in 
the  time  domain: 


_  3/2  „  d 

D.  (x)  =  - 


■1X2 


dt 


■1/2 


cx(t)  ] .  r  r(±y 


aj  I  ‘  f  x(D  (t-i)  ‘ir 


T{  l/a}  J  *  J  x  (  ^ 


♦  D*'2(x(0)) 


(37) 


Likewise  th«  one-half  derivative  can  be  obtained  from  either  x  or 
x,  and  it  is  defined  as: 

D;'*(x)  ■  *  [ru~>  ]-’  /  (M-lr 

+  D"2(x(0))  (38) 

The  initial  conditions  on  both  fractional  derivatives  are 
identically  zero  according  to  Bagley  (12).  The  other  three  time- 
dependent  quantities  in  the  equation: 

t 

x  =  J"  x  (T)dT  +  x(0)  (39) 

c 

t 

x  =  J  x(r)dT  +  x(0)  (40) 

c 

and  f  (t)  ,  an  arbitrary  time-dependent  forcing  function  which  is 
produced  in  a  separate  circuit  on  the  analog  computer. 
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In  addition  ail  terms  in  Eq  (35)  can  be  written  using  x  : 

( (t-  t)  _1/i  r 


x  =  m  1|^  f(t)  -  fll[  r(ixsj  J-1  J  x 
"  C2  [  J  x  <T)dT  +  *(0) j 


+  c3  ^  ru/z)  ]  V  [ .  Jx(T)dT+  ,(0)](t- 


,  -1/2, 
t)  d  T 


1  r  S. 

J  [  J  * 

5  to 


-  k  J  I  J'  x(T)dT  +  x(0)  I  dr 

o  L  o 


] 


(41) 


This  form  of  the  equation  motivates  modeling  this  second-order 
system  on  the  analog  computer.  That  is  one  operational 
amplifier  acts  as  x  and  is  a  summer  into  which  all  the  other  terms  serve 
as  input.  The  output  of  x  is  then  connected  to  an  integrator,  which 
produces  as  its  output,  x.  This  output  is  then  connected  to  another 
integrator,  and  its  output  is  x.  Use  of  a  half -order  fractional 

3/2  1/2 

differentiators  can  produce  the  (x)  and  Dt  (x)  as  previously 
explained.  Thu3,  all  the  terms  necessary  to  produce  x  are  available  and 
the  output  of  the  system  (x)  can  be  directed  to  a  display  or  recorder 
for  observation.  A  schematic  representation  of  this  circuit  is  found  in 
Figure  22. 

It  should  be  noted  the  operational  amplifiers  on  the  model  TB48 
are  inverting  amplifiers.  That  is  the  output  of  the  amplifier  has  the 
opposite  sign  of  the  input.  For  example,  when  x  is  put  through  an 
amplifier  configured  as  an  integrator,  the  output  is  -  x  instead  of  x. 

3/2 

To  obtain  the  correct  signs,  both  the  D  (x)  and  x  signals  must  be 
passed  through  inverters  to  obtain  the  correct  sign. 
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Figure  22.  General  Analog  Computer  Configuration  rOr 
a  Second-Order  System 
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Each  operational  amplifier  on  the  TR48  has  four  input  jacks;  two 
are  unity-gain  jacks  and  the  other  two  amplify  the  input  signal  by  a 
factor  of  ten.  These  factor  of  ten  jacks  will  hereafter  be  denoted  xlO 
To  provide  the  correct  coefficients  (c^  c2,  and  k)  the  output 
signal  for  each  derivative  term  must  be  connected  to  a  potentiometer 
(Figure  23) . 


Fiqure  23.  Potentiometer  Schematic 
and  Circuit  Symbol 
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Each  potentiometer  is  adjustable  in  i-'ioooth  increments  to 
provide  a  specific  fraction  of  the  amplitude  of  the  input  signal.  For 
example,  if  c  were  0.2500,  the  output  of  the  three-halves  derivative 
would  have  to  be  connected  to  a  potentiometer  set  at  0.2500.  However,  if 
c  were  1.2500,  the  three-halves  derivative  would  have  to  be  connected 
to  a  xlO  input  jack  on  the  x  summer.  Specific  details  for  establishing 
the  proper  signs  and  magnitudes  of  the  coefficients  are  presented  in 
Appendix  F. 

Open-  and  Closed  Loop  Simulation 

The  first  task  was  to  determine  a  set  of  coefficients  that  would 
result  in  frequencies  within  the  operational  band  of  the  half-order 
circuits  on  hand.  At  the  time  the  coefficients  were  developed,  the 
circuit  designed  by  the  continued  fraction  method  was  the  one  under 
consideration.  This  circuit  had  a  frequency  range  spanning  0.05  to 
0.40  Hz.  This  range  was  then  evaluated  to  determine  a  set  of  open- 
and  closed-loop  coefficients  compatible  with  the  circuit’s  performance. 

The  equation  developed  was: 

•  •  1/2 

x  +  2Dt  (x)  +  x  =  u(t)  (42) 

where  u(t)  is  the  control  force  (zero  for  the  open-loop  equation). 

It  should  be  noted  that  this  criteria  is  consistent  with  the 
Bagley-Torvik  model  for  a  viscoelastically-damped  structure  which  has  a 
stress-strain  relationship  modeled  by  the  addition  of  a  half-order  term 
(13:137) . 
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Using  the  linear  quadratic  feedback  theory  developed  in  Capt 
Walker's  thesis  (6),  the  optimal  feedback  for  this  system  was  determined 
to  be: 

u ( t )  =  1.6324  D*/2(x)  -  2.891  *  ♦  3.1556  D*'2(x)  -  0.4142x  (43) 

resulting  in  the  closed  loop  equation: 

x  -  1.6324  D2X2(x)  +  2.891  x  -  1.1556  D^2(x)  +  1.4142x  =  0  (44) 

Figure  24  illustrates  a  schematic  for  the  open-loop  equation,  and  Figure 
25  depicts  the  schematic  for  the  closed-loop  equation. 

The  initial  condition  on  x  for  both  the  open-  and  closed-loop 
cases  was  originally  set  at  2.5,  volts  and  in  the  final  stages  of  the 
experiment,  it  was  set  at  2.0  volts.  Stepping  down  the  -10  volt  reference 
voltage  on  the  attenuator  module,  through  a  potentiometer,  then  into  the 
IC  jack  on  the  x  operational  amplifier  established  this  initial 
condition.  This  voltage  acts  as  the  initial  disturbance  for  studying  the 
response  of  the  system.  Again,  referencing  Figures  24  and  25,  both 
half-order  circuits  are  isolated  from  the  initial  conditions  by- 
comparators.  When  the  'OP*  button  is  pushed  to  start  the  simulation  the 
comparators  connect  the  input  to  the  one-half  and  three-halves-order 
circuits.  This  was  accomplished  to  conform  to  the  necessary  initial 
>  conditions  established  by  Bagley  for  all  fractional  order  derivatives 

(12).  The  analog  computer  was  operated  according  to  the  detailed 


i 
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description  discussed  in  Appendix  F,  and  the  output  is  recorded  on  the 
Hewlett-Packard  recorder  (model  7090A) .  Figure  26  is  a  sample  output; 
it  illustrates  the  x  response  for  both  the  open-  and  closed-loop 
systems . 
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Figure  25.  Closed-Loop  System  with  Comparator  Isolation 
o/Of(x)  Input 


Two  theoretical  methods  for  predicting  the  open-  and 
closed-loop  responses  were  developed  in  parallel  with  the  experimental 
effort.  One  method  involved  Laplace  transforms,  and  it  is  detailed  in 
Appendix  A.  The  other  uses  a  Miitag-Lef f ler  expansion  in  the  time 
domain,  and  it  is  detailed  in  Appendix  C.  Both  methods  produced  exactly 
the  same  response  for  the  open-  and  closed-loop  cases.  Therefore,  only 
the  Laplace  transform  response  will  be  shown  for  comparison  (Figure  27 
for  the  open-loop  and  Figure  28  for  the  closed-loop  response) .  It  is 
easy  to  recognize  the  initial  discrepancies  between  the  experimental 
results  and  the  theoretical  predictions.  In  the  open-loop  case  the 
analytic  methods  predict  a  very  heavily  damped  response  that  never 
crosses  the  horizontal  axis  and  approaches  zero  from  above.  The 
experimental  results  crossed  the  horizontal  axis,  approached  zero  from 
below  and  decayed  rapidly.  The  same  type  of  disparity  was  observed  on 
the  closed-loop  system.  Although  feedback  made  the  system  more 
responsive,  it  became  necessary  to  resolve  the  discrepancy  between 
theory  and  experiment. 

Total-Cycle  Simulation 

The  previous  con' iguration  modeled  the  open-  and  closed-loop  cases 
as  initial-value  problems.  This  procedure  starts  with  the  system  at 
rest.  By  using  a  function  switch  (described  in  Appendix  F) ,  a  step 
voltage  was  applied  to  the  circuit  (Figure  29). 
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Figure  28.  Theoretical  Closed-Loop  Response 
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Figure  29.  Total— Cycle  Analog  Configuration 


The  circuit  undergoes  a  transient  response,  eventually  coming 
to  rest  at  a  displacement  value  approximately  equal  in  magnitude  to  the 
size  of  the  step  function.  This  steady  system  state  is  nearly  the  same 
as  that  modeled  in  the  original  configuration  by  placing  an  initial 
condition  on  x.  The  one  difference  in  this  case  is  that  the 
fractional-order  derivatives  are  not  isolated  from  their  inputs. 

When  the  function  switch  was  toggled,  the  step  function  was  removed. This 
section  of  the  response  curve  approximates  the  response  modeled  by  the 
theoretical  predictions.  The  result  of  this  total-cycle  simulation  is 
shown  in  Figure  30.  Figure  31,  for  the  open-loop  case,  and  Figure  32, 
for  the  closed-loop  case,  represent  the  second  half  of  the  Total-Cycle 
Simulation.  The  experimental  response  compares  favorably  with  that 
modeled  by  the  initial-value  problems.  Section  V  (Analysis)  draws 
conclusions  concerning  the  modifications  required  to  the  ini tial -value 
problem  configuration  to  more  correctly  model  the  initial  conditions. 
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Final  Configuration  Performance  Evaluation 

The  final  configurations  for  both  open-  and  closed-loop  initial- 
value  problems  are  pictured  in  Figures  33  and  34,  respectively.  There 
are  two  changes  made  from  the  original  initial  value  configuration. 
First,  both  half-order  circuits  are  not  isolated,  but  allowed  to  charge 
their  capacitors  to  whatever  value  is  consistent  with  their  input. 
Secondly  the  x  summer  was  isolated  from  its  input  using  a  comparator 
(Appendix  F)  until  the  ‘OP*  button  was  pressed.  Final  results  for  both 
open-  and  closed-loop  cases  are  compared  with  the  theoretical 
predictions  in  Figures  35  and  36.  The  minor  discrepancies  are  discussed 
in  the  analysis  section. 
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Figure  34.  Final  Closed-Loop  Configuration 
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V.  Analysis  of  Results 


Introduction 

This  section  deals  only  with  anomalies  and  unexpected  results 
encountered  during  the  research.  Tasks  which  closely  matched 
theoretical  predictions  are  covered  in  the  conclusion  section;  it  is 
sufficient  tc  state  they  functioned  as  expected.  The  discussion  is 
divided  into  three  parts;  problems  in  the  circuit  validation,  problems 
with  simulation  results,  and  disadvantages  associated  with  using  the 
analytical  tools.  In  each  area,  the  problems  are  defined,  possible 
causes  listed,  and  solutions  presented. 

Circuit  Validation 

The  voltage  level  output  from  the  Tektronix  function  generator  was 
not  constant  over  the  frequency  range  tested.  At  low  frequencies,  (0.05 
Hz)  the  output  voltage  was  several  volts  lower  than  at  higher 
frequencies.  The  function  generator  is  designed  to  terminate  into  a  500 
impedance  load.  The  Oldham-Zoski  circuit  impedance  is  much  higher  than 
this,  especially  at  low  frequencies  where  the  majority  of  the  current 
passes  through  the  resistors.  By  placing  a  500  resistor  in  parallel 
with  the  circuit,  the  impedance  remained  close  to  500  .  The  variation 
of  circuit  impedance  with  the  frequency  at  which  the  circuit  was  tested 
was  not  observed  to  adversely  affect  the  performance  of  the  test 
configuration.  Slight  diiferences  were  noted  in  the  start-up  transient 
response  of  the  input  and  output  signals  as  a  function  of  frequency. 
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These  transient  characteristics  are  repeatable  for  a  specific  circuit 
and  frequency.  In  addition,  the  amplitude  of  the  input  signal  remained 
between  8.00  and  8.10  volts  for  all  frequencies  tested  (Appendix  G 
contains  the  circuit  verification  data).  These  slight  variations  were 
judged  to  be  acceptable. 

When  using  the  oscilloscope  and  counter/timer  combination  to 
determine  magnitude  and  phase  angle,  the  results  were  not 
repeatable  to  the  accuracy  desired  (±1  dB/decade  gain  slope  and  ±5 
degrees  phase  shift).  The  slope  of  the  Bode  magnitude  plot  would  vary 
from  9.1  to  11.25  dB/decade,  and  the  phase  angle  varied  from  37  to 
43  degrees  as  a  function  of  frequency.  The  oscilloscope  was  not 
capable  of  tracking  the  output  of  the  fractional  order  circuit  using 
just  one  magnitude  range.  Constant  adjustment  was  required.  It  became 
necessary  to  check  the  adjustments  twice  at  each  frequency  to  ensure 
that  the  controls  had  not  been  bumped  when  changing  from  one  channel  to 
the  other.  The  trigger  control  on  the  counter/timer  which  is  tuned  to 
detect  a  signal  crossing  through  the  zero  voltage  level  was  extremely 
sensitive.  Because  of  the  need  for  frequent  adjustments,  the  confidence 
in  these  results  was  low. 

It  became  necessary  to  attempt  verification  of  the  circuit 
on  the  analog  computer.  The  three-halves  derivative  functioned  correctly 
on  the  analog,  when  implemented.  The  phase  was  approximately  43 
degrees,  and  the  magnitude  could  be  adjusted  correctly.  However, 
running  the  output  of  the  three-halves  derivative  through  an  integrator  to 
produce  the  one-half  derivative  produced  the  sinusoid-ramp  combination 
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illustrated  in  Figure  17.  The  peak-to-peak  amplitude  of  the  sine  wave 
had  the  correct  magnitude,  but  the  underlying  ramp  function  indicated  a 
bias  voltage  at  the  input  of  the  operational  amplifier. 

The  source  of  the  bias  could  not  be  readily  identified.  With 
visual  observation  of  the  phase  and  magnitude,  it  was  impossible  to 
determine  small  biases  on  the  output  signal  or  to  observe  the  transient 
adjustment  of  the  circuit  to  an  input  signal.  To  further  complicate 
matters,  the  original  method  of  connecting  the  fractional-order  circuit 
to  the  operational  amplifier  was  accomplished  by  inserting  a  banana  plug 
directly  into  the  summing  junction  of  the  operational  amplifier.  This 
procedure  avoided  adding  resistance  to  the  circuit  by  connecting  into 
one  of  the  input  resistors  of  the  operational  amplifier.  When  the 
signal  passed  through  the  operational  amplifier,  it  produced  the 
extremely  noisy  output  signal  illustrated  in  Figure  37. 

Many  possible  causes  were  identified  with  the  original  circuit 
configuration.  First,  there  was  concern  with  the  sirv  potential  sources 
of  stray  voltages  and  inductances.  The  operational  amplifiers  were  not 
commonly  grounded;  this  situation  enhanced  possible  ground  potential 
biases  between  components.  The  long  patch  cords  connecting  the  circuit 
board  to  the  analog  computer  passed  in  close  proximity  to  the  analog 
computer's  power  supply.  The  long  leads  on  the  discrete  resistors  and 
capacitors  were  not  trimmed,  resulting  in  a  circuit  similar  to  that 
depicted  in  Figure  38.  These  items  were  corrected  by  grounding  the 
circuit  board  and  operational  amplifiers  together,  trimming  component 
leads  so  that  the  components  were  flush-mounted  to  the  circuit-board, 
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Figure  38.  Long  Leads  on  Discrete  Components 

and  using  coax  cables  to  connect  the  circuit-board  to  the  analog 
computer.  Reducing  the  final  resistance  value  in  the  half-order  circuit 
by  10KCI  and  patching  its  output  into  a  X10  input  jack  (10K^  resistance) 
on  the  operational  amplifier,  instead  of  the  summing  junction,  eliminated 
the  noisy  signal . 

The  most  significant  change  was  the  use  of  the  Hewlett-Packard 
analog-digital  recorder  (model  7090A)  It  eliminated  problems  associated 
with  the  start-up  transient.  Appendix  G  contains  the  plots  of  the 
circuit  validation;  the  start-up  transient  is  seen  to  be  smooth  with  no 
discontinuities.  The  analog-digital  recorder  greatly  simplified  data 
collection  and  analysis;  set-ups  became  routinely  simple,  and  digitized 
values  of  interest  on  the  output  plot  could  be  annotated  directly  on  the 
plot  through  the  use  of  a  built-in  routine. 
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The  method  of  supplying  the  one-half  derivative  was  altered  to 
eliminate  the  ramp.  A  second  half-order  circuit  was  constructed  with 
capacitors  having  twice  the  value  of  the  capacitors  used  in  the  first 
circuit  and  resistors  having  half  the  value  of  those  used  in  the  original 
circuit  (the  time  constant  remaining  the  same) . 

The  performance  comparisons  of  the  circuits  produced  an  unexpected 
result.  When  second-order  effects  are  considered,  a  capacitor  is 
modeled  as  a  lossless  capacitance  in  parallel  with  a  resistance  (19). 

When  combined  into  a  typical  Oldham-Zoski  cell,  the  second-order  model 
then  has  a  capacitance  in  parallel  with  two  parallel  resistances.  One 
would  predict  the  circuit  with  the  lowest  capacitance  to  be  the  circuit 
with  the  best  phase  characteristics.  Comparison  of  Figures  7  and  9  show 
this  is  indeed  the  case.  When  designing  future  circuits  an  effort 
should  be  made  to  make  the  capacitor  values  as  small  as  feasible.  The 
use  of  an  impedance  bridge  to  measure  the  'quality*  of  the  capacitances 
is  also  recommended. 

The  circuit  produced  a  phase  angle  of  43  degrees  when  connected 
to  the  analog  computer,  but  only  37-38  degrees  when  measured  by 
oscilloscope  and  analog-digital  recorder.  Examination  of  the  input 
circuitry  on  the  analog-digital  recorder  showed  a  lMft  resistor  connected 
in  series  with  the  input.  To  correct  this  situation  the  circuit  was 
connected  to  ground,  instead  of  connecting  the  half-order  circuit 
directly  to  the  input  of  the  recorder.  The  recorder  was  then  put  in 
parallel  with  the  last  resistor  of  the  half-order  circuit.  This 
configuration  produced  a  phase  angle  of  43  degrees. 


Finally,  the  ramp  in  the  one-half  derivative  that  resulted  from 
integrating  the  three-halves  derivatives  requires  analysis.  With  the 
Hewlett-Packard  analog-digital  recorder,  small  amplitude  biases  were 
observed  in  the  output  of  both  half-order  circuits.  The  input  from  the 
function  generator  was  also  observed  to  be  biased.  It  appeared  this  was 
the  reason  for  the  anomaly.  A  closer  evaluation  proved  differently. 

For  example,  in  the  validation  plot  (Appendix  G)  for  Circuit  1  at 
f  =  0.1451Hz  (5;  to  s  i) f  the  ratio  of  the  input  amplitude  (peak-to-peak) 
to  the  output  amplitude  is  63.32.  However,  the  ratio  of  total  offset 
is: 
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The  bias  observed  was  actually  ouch  less  than  that  predicted  by 
this  first-order  method.  It  is  noted  that  Circuit  1  is  biased  positively 
while  the  one-half  derivative  exhibits  a  negative  bias.  The  operational 
amplifiers  involved  are  inverters;  therefore,  the  bias  does  have  the 
correct  sign. 

Does  this  bias  affect  the  simulation,  and  is  it  an  error  that 
increases  with  time?  The  output  of  the  one-half  and  three-halves 
derivatives  is  connected  to  the  summer  and  results  in  an  amplitude  error 
of  the  opposite  sign  of  the  particular  fractional  derivative;  that  is 

3/2  X/2 

x  s  -d^  (x)  -  Dt  (x)  +  .  .  .  .  If  the  one-half  derivative  is  larger 

than  nominally  predicted,  and  it  should  be,  then  x  is  slightly  smaller 
than  it  should  be.  Thus  x  is  smaller  then  it  would  be  if  the  bias  was 
not  there,  and  the  same  argument  applies  to  x.  When  x  is  differentiated 

X/2 

by  the  fractional-order  circuit  the  result  is  a  slighty  smaller  (x) 
input,  compared  to  the  case  where  the  bias  was  not  present.  In  other 
words,  with  feedback,  the  error  is  stable  and  does  not  grow  with  time. 

Upon  completion  of  the  second  fractional-  order  circuit  design,  the 
circuit  was  connected  to  the  harmonic  oscillator  for  an  amplitude  check 
with  oi  =  l  radian/sec.  The  harmonic  oscillator’s  output  amplitude  was 
not  stable,  and  decayed  with  time.  The  first  thought  was  to  remove  the 
fractional-order  circuit  to  see  if  it  might  be  the  cause.  With  this 
done,  the  amplitude  continued  to  decay.  Three  new  operational 
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amplifiers  were  selected,  the  circuit  was  rewired,  and  the  problem 
persisted,  indicating  it  was  not  an  operational  amplifier  problem.  A 
second  analog  computer  was  placed  into  service,  thinking  the  power 
supply  was  bad  on  the  first  unit,  but  the  decay  still  continued.  It  was 
then  noticed  that  everytime  someone  used  a  computer-printer  in  the  room, 
the  voltage  would  experience  a  step-function  loss  of  amplitude.  Waiting 
until  the  demand  decreased  on  the  electrical  supply  circuit  before 
operating  the  analog  computer  solved  the  problem.  The  whole  analog 
computer  functioned  better;  the  potentiometers  were  not  as  noisy  and 
easier  to  set.  The  operating  schedule  was  shifted  to  occur  late  at 
night  when  the  electrical  supply  was  more  stable. 

Open-  and  Closed-Loop  Simulation  Problems 

Figure  26  shows  the  initial  open-  and  closed-loop  responses  obtained 
experimentally.  Figures  27  and  28  depict  the  results  predicted  for  the 
simulation  by  the  Laplace  transform  method.  The  frequency,  amplitude, 
and  axis  crossings  vary  significantly  from  that  observed  experimentally. 

A  second  Laplace  transform  prediction  was  made  (Figure  39) ,  deleting  the 
s  terms  when  transforming  the  half-order  derivatives  (see  Appendix 
A).  This  response  was  a  almost  entirely  due  to  the  residue,  there  was 
very  little  contribution  from  the  integral  along  the  branch  cut.  The 
result  was  in  close  agreement  with  the  experimental  data,  however,  it 
did  not  agree  with  the  initial-value  problem  work  reported  by  Bagley 

-  i/2 

(12)  which  required  the  s  terms. 
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RESPONSE  (VOLTS) 


Figure  39.  Theoretical  Response  Without  the  s-l/2 
Initial-Value  terms 


The  first  thought  was  that  the  circuit  and  experimental  methods 

_  £/2 

were  incorrect.  The  second  thought  was  that  the  s  terms  in  the 

Laplace  transform  method  were  invalid.  After  looking  at  several 
possibilities ,  it  appeared  that  the  problem  centered  on  the  method  of 
applying  the  initial  conditio"*  and  an  understanding  of  the  definition  of 
a  zero  initial  condition  on  the  half-order  circuits. 

Two  specific  applications  of  initial  conditions  were  in  question. 
During  the  gain  evaluation  phase  a  decision  was  made  to  isolate 
the  one-half  derivative’s  input  until  the  'OP'  button  was  pressed. 

Figure  40  shows  the  response  of  the  one-half  derivative  when  not 
isolated.  It  appears  to  require  four  or  five  cycles  to  adjust  to  the 
correct  amplitude.  Comparison  of  this  response  with  that  of  Figure  41, 
in  which  the  one-half  derivative  was  isolated  with  a  comparator,  led  us 
to  believe  that  a  zero  initial  condition  on  the  half-order  circuits 
meant  isolation.  When  isolated,  the  one-half  derivative  attains  the 
correct  amplitude  in  one  cycle.  For  this  reason  the  decision  was  made 
to  isolate  both  the  three-halves  and  one-half  derivatives.  Secondly, 
when  a  system  is  at  rest  at  some  displacement  value,  it  is  expected  that 
the  acceleration  would  be  zero.  In  the  simulation,  the  x  initial 
condition  of  2.0  volts  appeared  a3  -2.0  volts  at  the  summer,  indicating 
the  system  was  not  at  rest. 

Two  te3t  cases  were  devised  to  address  the  problem.  In  both 
cases,  the  initial  displacement  was  set  to  zero.  In  the  first  case,  a 
harmonic  oscillator  with  to  =  1.585  rad/sec  served  as  the  input  to  the  x 
summer  as  a  forcing  function  for  the  circuit  x  +  x  =  0.  The  one-half 


derivative  was  connected  to  the  displacement  (x)  output.  The  response 

of  the  one-half  derivative  and  displacement  were  recorded. 

1/2 

If  a  45  degree  phase  lag  and  an  amplitude  of  (1.585)  was  observed 
for  the  one-half  derivative,  the  circuit  was  working  correctly.  After  one 
complete  cycle,  this  was  the  case  (Figure  42).  This  test  verified  the 
circuit  and  methods  used;  therefore,  the  initial  required  further 
investigation.  The  comparators  were  removed  from  the  half  derivatives 
and  x(0)  was  set  equal  to  zero.  The  "Total-Cycle  simulation",  as 
described  in  the  Section  IV,  was  then  run.  Figures  30  and  31  are 
examples  of  this  procedure  for  the  open-loop  system  and  Figure  32  is  an 
example  for  the  closed-loop.  On  cne  of  the  runs  the  "HD"  (hold)  key  was 
used  to  hold  the  values  of  the  displacement  (x)  operational  amplifier  at 
the  upper  stationary  position.  The  summer  value  displayed  on  the 
Digital  Volt  Meter  (DVM)  at  this  point  was  0.0000  volts.  With  the 
circuit  still  on  hold,  the  function  switch  was  then  opened,  the  value  of 
x  changed  immediately  to  the  negative  of  the  value  of  the  step  function. 
The  ‘OP"  button  was  pressed  and  the  circuit  duplicated  the  second  half 
of  the  total-cycle  simulation.  Armed  with  this  information  x  was 
isolated  from  all  inputs  until  pressing  the  "OP"  button  closed  the 
comparator  contacts.  Neither  of  the  half-order  circuits  was  isolated. 

The  circuit  was  reconfigured  so  the  initial  condition  was  again  applied 
through  the  IC  jack  in  the  x  operational  amplifier.  Figures  35  and  36 
depict  the  final  experimental/analytic  comparisons  for  both  the  open- 
and  closed-loop  simulations. 
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Figure  40.  Half  Derivative  Response  with  the  Input 
Not  Isolated  at  t=0 


Figure  42.  Response  to  a  Sinusoidal  Forcing  Function 
with  (m2  =3.5 


It  appears  that  to  function  properly  in  an  initial-value 
problem  simulation  the  capacitors  in  the  fractional-order  circuits  must 
be  charged  to  a  voltage  consistent  with  the  initial  value  of  their  input 
signal.  If  the  initial  value  of  the  input  to  the  half-order  circuit  is 
zero,  zero  initial  voltage  on  the  half-order  circuit  is  consistent.  If 
the  initial  value  of  the  signal  input  to  the  fractional -order  circuit  is 
other  than  zero,  then  the  initial  voltage  of  the  circuit  must  assume  a 
voltage  consistent  with  that  input  voltage. 

There  was  a  small,  apparently  constant,  displacement  in  the 

experimental  value  when  compared  with  the  analytical  results.  If  the 

absolute  value  of  the  displacement  varied  as  the  ratio  of  c  to  k  in  the 

equations,  it  would  indicate  a  relationship  between  the  displacement  and 

one-half  derivative.  However,  the  magnitude  of  the  displacement  does 

not  seem  to  vary  from  the  open-loop  to  the  closed-loop  situation;  it 

1  /2 

only  changes  signs.  When  not  isolated,  the  (x)  initial  value  is 
0.143  volts  in  the  open-loop  case  and  0.135  volts  in  the  closed-loop 
configuration.  The  net  displacement  between  the  experimental  and 
analytical  results  is  approximately  0.07  volts.  This  value  is  nearly 
the  same  magnitude  as  the  circuit  bias,  but  if  the  bias  is  the  cause, 
the  value  of  the  displacement  should  change  by  a  factor  of  two  due  to 
the  change  in  c  .  However,  the  0Z  circuit  is  only  an  approximation  to 

1X2 

an  s  transfer  function  over  the  frequency  range  of  its  design. 
Therefore,  the  long  term  accuracy  of  the  circuit  is  dependent  on  its 
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lower  frequency  limit.  Consequently,  the  results  observed  are  likely  to 


be  the  best  possible  for  a  circuit  with  this  set  of  design  parameters. 
These  displacements,  in  both  cases,  make  the  circuit  more  responsive  and 
approach  zero  faster  than  predicted .This  anomaly  will  be  left  as  an  open 
issue  for  follow-on  research. 

Theoretical  Tools 

Mittag-Lef fler  Expansion  (Appendices  C  and  D) 

The  program  MITLEFR.DAT,  an  executable  file  in  Matrix^, 

expands  the  roots  of  the  characteristic  equation  in  an  infinite  series 

to  provide  the  response  of  the  open-and  closed-loop  systems  in  the  time 

domain.  For  this  case  the  solution  was  obtained  for  a  second-order 

1/2 

system  involving  multiples  of  the  half  derivative  [x.D^  (x)  ,  x, 

3/2 

(x)  ].  However,  it  can  be  used  to  expand  a  second-order  equation 

made  up  of  any  1/nth  fractional-order  terms  by  entering  (N  =  2xn)  when 

prompted.  Input  parameters  are  upper  and  lower  time  limits  for 

calculation,  actual  response  start  time,  coefficients  of  the 

1  ^2 

differential  equation,  and  initial  conditions  for  xf  *,Dt  (x) ,  and 

3/2 

(x) .  The  method  is  effective  for  approximately  15  seconds  from  the 
time  the  response  starts.  For  extended  lengths  of  time,  the  number 
of  terms  required  to  obtain  convergence  of  the  expansion  becomes 
excessive.  For  example,  it  required  2  hours  of  run  time  to  evaluate  100 
points  using  200  terms  in  the  series  for  the  open-loop  system.  It 
required  nearly  the  same  amount  of  time  to  evaluats  50  points  using  400 
terms  in  the  closed-loop  case 


(  86  ) 


Laplace  Transform  (Appendices  A  and  B) 


The  program  RESPONSE  is  written  in  Fortran  and  performs  a 
contour  integration  in  the  complex  plane  to  calculate  the  time  response 
(20:818-824).  The  current  method  requires  manual  calculation  of  the 
residues.  In  addition,  asymptotic  approximations  to  the  integrand  must 
be  made  to  simplify  the  inversion  process.  Figure  43  shows  a  comparison 
of  five  different  selections  for  the  asymptotic  limits.  Three  of  the 
five  trials  required  twenty  minutes  of  computational  time;  the  other  two 
required  two  hours. 

The  results  of  the  five  trials  are  remarkably  close,  yet  the 
processing  time  can  be  increased  unnecessarily  by  an  improper  selection 
of  limits.  Each  time  an  analysis  is  accomplished  for  a  different  set  of 
coefficients,  the  equations  for  four  derivatives  of  three  functions  must 
be  inserted  into  the  proper  function  routine  in  the  program.  The 
supporting  function  routines  must  be  recompiled  and  then  linked  with  the 
main  program  for  each  specific  set  of  equation  coefficients  (m,c  , 
c  ,k) .  The  Laplace  transform  method  is  not  limited  with  respect  to 

n 

time.  If  it  is  necessary  to  characterize  the  entire  response  of  a 
specific  system,  and  that  system’s  response  lasts  longer  than  15 
seconds,  the  Laplace  method  is  viable.  However,  the  Laplace  transform 
method  is  limited  with  respect  to  responses  which  have  Laplace 
transforms . 

Additional  Commments  -  Deciding  Which  Method  to  Choose 

Figure  44  shows  a  comparison  of  the  two  methods;  they  are  identical. 
For  initial  characterization  of  a  system  the  Mittag-Lef f ler  should 
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be  used  since  it’s  easy  to  implement.  For  extensive  time-domain 
analysis  the  Laplace  transform  method  is  preferred. 

Why  use  the  analytic  tools,  if  the  analog  computer  i3  available’ 

It  was  the  Laplace  transform  method  which  identified  the  problem  with 
invalid  initial  conditions.  Figure  45  illustrates  a  preliminary  analog 
simulation  of  the  closed  loop-system.  The  sign  on  the  three-halves 
derivative  was  incorrect.  Had  it  not  been  for  the  analytical 
prediction,  this  situation  would  not  have  been  identified.  In  the 
future,  if  someone  develops  a  generalized  Laplace  transform  method  that 
solves  the  problem  in  terms  of  the  coefficients  (m,k . c  ),  receives 

n 

them  as  input  parameters,  and  then  calculates  the  residues,  the  solution 
for  a  given  system  could  be  available  in  twenty  minutes.  This  tool  is 
probably  quicker  than  the  analog  simulation  technique.  However,  both 
analog  and  theoretical  methods  have  special  tasks  that  they’re  more 
efficient  at.  It’s  best  to  improve  both  methods. 


RESPONSE  (VOLTS) 


2.50 
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Figure  45.  Circuit  Response  with  Improper  Coefficient 


VI  Conclusions  and  Recommendations 


I 


I 
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Summary 

This  thesis  demonstrated  that  analog  simulation  of  differential 
equations  involving  fractional-order  derivatives  can  be  performed. 
Fractional-order  feedback  in  linear  systems  is  feasible,  and  the  response 
is  predictable  with  computer  methods.  This  feedback  is  stable  and 
improves  system  performance.  The  methods  developed  are  directly 
applicable  to  systems  whose  behavior  can  best  be  described  by  fractional 
derivatives.  These  methods  could  provide  additional  degrees  of  freedom 
for  control  engineers  to  improve  system  performance  (12).  With 
the  additional  fractional-order  time  derivatives  available  for  feedback, 
response  can  be  tailored  to  precisely  fit  the  needs  of  the  application 
without  a  loss  of  system  stability. 

The  Oldham-Zoski  circuit  provides  the  l^nth  fractional  derivative 
of  an  input  signal.  Previous  research  indicated  the  circuit  functioned 
for  constant  and  ramp  input  signals,  for  any  fractional  integral  cr 
derivative.  This  investigation  demonstrated  that  this  circuit  also 
performs  with  a  sinusoidal  input.  Therefore,  it  is  anticipated  that 
general  fractional-order  circuits  can  be  designed  and  fabricated  to 
cover  the  frequency  range  and  functional  type  required  by  a  specific 
application. 


» 
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The  circuits  designed  and  fabricated  can  not  only  function  as 
tools  to  simulate  processes  involving  fractional  derivatives,  but  they 
can  also  be  used  to  modify  traditional  sensors  for  active  control 
systems . 

Two  analytic  tools  were  developed  to  predict  the  experimental 
response.  Fully  functional  computer  programs  were  developed  for  both 
methods.  Even  though  the  two  methods  were  totally  unrelated,  the 
results  generated  were  indistinguishable.  A  high  degree  of  confidence 
is  merited  for  these  analytic  tools  as  they  not  only  described  the 
response,  but  they  also  predicted  it. 

Areas  for  Further  Study  and  Development 

1)  Perform  simulation  and  feedback  including  higher  orders  of 
fractional  derivatives  as  well  as  the  one-half  and  three-halves. 

Compare  the  response  to  the  results  presented  here.  Investigate  if 
the  added  degrees  of  freedom  produce  better  system  response. 

2)  Electrical  circuit  development.  Investigate  designing 

and  fabricating  the  fractional -order  circuits  using  current  IC  methods. 
Determine  if  the  circuit  bias  is  inherent  in  the  design  or  can  be 
removed.  If  it  can  be  compensated  for,  a  single  fractional-order 
micro-circuit  could  supply  the  l^nth  and  n+i/nth  derivative.  In  addition, 
define  interface  requirements  for  the  circuits,  that  is  power  consumption, 
impedance  matching,  etc.  If  added  degrees  of  freedom  are  beneficial, 
determine  the  number  that  can  be  designed  into  one  integrated  circuit. 


3)  Investigate  the  steady-state  offset  of  the  experimental  values. 
Determine  if  a  designing  a  circuit  for  operating  at  lower  frequencies 
produces  an  improved  steady-state  response  or  if  the  steady-state  offset 
observed  was  a  function  of  the  initial  charge  on  the  capacitors.  One 
possibility  of  testing  for  this  is  to  run  a  simulation  with  a  higher 
initial  charge  on  the  capacitors,  and  then  observe  the  results. 

4)  Continue  development  of  the  Laplace  inversion  program. 

Automate  the  calculation  of  residues.  Solve  in  terms  of  general 
coefficients  (k.m.c^,  ,..,cn);  use  the  general  coefficients  as  input  to 
the  function  routines. 

5)  Combine  the  circuit  with  an  active  control  system  and  logic  to 
control  the  response  of  an  actual  structure.  One  possibility  -  see  if 
the  AFIT/ENG  (Department  of  Electrical  and  Computer  Engineering)  robotics 
laboratory  is  interested. 

6)  Investigate  the  practical  limits  of  frequency  range  available. 
Determine  if  a  circuit  designed  to  operate  over  two  decades  of  frequency 
functions  as  effectively  over  those  two  decades  as  one  designed  for  six. 

7)  Within  the  design  frequency  range  of  the  circuit  the  transfer 
function’s  slope  is  10  dB/decade  with  a  45  degree  phase  shift. 

Determine  the  operating  characteristics  of  the  circuit  outside  its 
design  range. 

8)  Without  scaling  on  the  analog  computer  the  voltage  limits 

restrict  the  magnitudes  of  voltage  that  can  be  used  at  high  frequencies. 

* 

Research  the  necessary  modifications  to  the  circuit  for  t  operating 
* 

times,  where  t  is  a  fraction  of  the  actual  time. 
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Equipment  Improvements 

1)  A  Hewlett-Packard,  or  similar  PC,  with  an  IEEE-488-bus 
compatible  interface  card  for  connecting  to  the  Hewlett-Packard 
analog-digital  recorder  (model  7090A) .  Data  from  the  recorder  could  be 
directly  transferred  to  the  PC  for  analysis.  Likewise  theoretical  plots 
could  be  done  on  the  recorder.  The  Hewlett-Packard  recorder  (model 
7090A)  has  additional  capabilities  when  it  is  interfaced  with  a  PC  than 
it  has  in  its  stand-alone  mode. 

2)  If  it  is  determined  to  still  be  more  desirable  to  use  discrete 
components,  obtain  a  supply  of  IX  components. 

3)  Investigate  the  use  of  the  AFIT  Department  of  Electrical 
and  Computer  Engineering  (AFIT/EMG)  hybrid  analog  computer  (SIMSTAR) . 

It  is  the  latest  model  available,  can  be  programmed  through  a  terminal 
using  Fortran  -  as  opposed  to  using  patchcords,  and  has  an  interface  for 
external  circuitry. 
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Laplace  Transform  Method  for  Theoretically  Predicting  Time  Response 


A  general  second-order  equation  for  a  viscoelastical ly  damped 
system  can  be  written: 

"  3/2  1/2 

mx  +  (x)  +  c  £ '  +  c  J)  t  lx)  +  kx  =  f(t)  (Al) 

where  x,  x,  and  x  represent  acceleration,  velocity,  and  position. 

3/2  1/2 

(x)  and  Dt  (x)  are  the  three-halves  and  one-half  fractional 
derivatives.  These  fractional  derivatives  aid  in  providing  simple, 
causal  models  of  the  behavior  of  analytic  systems  (2:125-126).  (  mx 

is  an  inertial  term,  c2x  is  viscous  damping,  and  kx  is  stiffness. 

3/2 

Therefore,  (x)  can  be  defined  as  a  viscoinertial  damping  term,  and 
1/2 

Dt  (x)  as  a  viscoelastic  damping  term.) 

Open-Loop  Example  Case 

For  an  example  case,  let  us  examine  a  model  of  a  typical 
viscoelastically  damped  system; 

c  =  c2=  f(t)  =  0  ,  m  =  k  =  1  ,  cg  =  2: 

x  +  2D*Xi(x)  +  x  =  0  (A2) 

The  initial  conditions  are: 

x (0)  =  A  and  x(0)  =  0  (A3) 


To  solve  the  differential  equation,  transfer  it  to  the  Laplace  domain. 
Integer-order  derivatives  transform  as  usual,  but  fractional-order 


derivatives  transform  slightly  differently.  The  mathematical  definition 
of  an  m/n  -order  derivative  contains  an  integral  expression: 


d‘x2(x)  =  ]  ~Af  -  d  t 

t  J  1/2 

C  T 


and  in  the  Laplace  domain  (12): 


£{  dJ/2(x)>  =  s'^s*  -  x (0)  1  =  s1^2-  s"1/2(0) 


and  the  Laplace  transform  of  Eq  (A2)  becomes: 


-  As  +  2 a"**  -2A S_1X2  +  x  =  0 


and  solving  for  IK(s) : 


(s  +  2s~1/Z) 
s2  +  2s1XZ  +  1 


To  obtain  the  response  of  the  system,  one  must  calculate  the 
inverse  Laplace  transform  of  )K(s)  (29:818-824,14:141-143): 

J'+i.OO 

£~X{  ^(s))  =  — —  f  i*((s)eat as  (A8) 

2TTv  J 

vOO 

Solving  Eq  (A8)  involves  integrating  around  a  closed  contour,  with 

1/2 

Eq  (A8)  as  one  of  the  contours  (Figure  46).  Again,  the  s  terms 

make  the  solution  technique  slightly  different  than  usual.  Mapping  the 

*  A  1/2 

denominator  of  Eq  (A7)  into  the  s-plane,  where  s  =  s  ,  permits  the 

*  4 

characteristic  equation  to  be  written  s  +  2s+  1.  The  contour  integral 

1/2 

can  then  either  be  performed  in  the  s  -plane,  or  remapped  back  into 


the  s-plane  (21:197-199) . 


(97) 


Figure  46.  Contour  Integral  Definition  for  the 
Laplace  Transform  Method 


(98) 


Pia  _1*  p2=  *0-5437-  P3P3=  0-7718  ±11.1151  (A9) 

These  are  the  poles,  and  are  then  mapped  back  into  the  s-plane  by  the 
transformation : 


s  =  o  +  iu 
s  =  u  +  iv 

"  i 

s  =  s 

2  2 

=  u  -  v  +  2iv 

Equate  the  real  and  imaginary  parts: 

2  2 

o  =  u  -  v  ;  w  =  2iv 


(A10) 


(All) 


The  poles  in  the  s  plane  are  now: 

pt=  1.0,  p2=  0.29561,  p3p"=  -0.64777  ±il. 72127  (A12) 

Normally,  poles  on  the  positive  real  axis  indicate  an  unstable  system. 
However,  in  this  case,  they  result  from  the  s-to-s  mapping.  These  poles 
originated  in  the  area  of  the  s-plane  which  did  not  map  onto  the 
principle  sheet  of  the  Riemann  surface  in  the  s  plane.  Their  effect  on 
the  system  comes  from  integrating  along  the  branch  cut  depicted  in 
Figure  46. 

The  mathematical  expression  for  the  residue  theorem  is  (20:818): 

“5ST  (l  +  J  +  J  +  I+  J+  J  )  Z  Residue  <v>  (A13) 

'•l  2  3  4  5  6  J  1  =  1 

The  Laplace  inversion  integral,  Eq  (A8) ,  is  represented  by 
integration  along  path  1.  Therefore: 

X(.)  )  =  +  2  Residue  <v>  (A14) 

'  2  3  3  4  5  v  =  i 


The  quantity  on  the  left  hand  side  of  Eq  (A14)  is  defined  by  Eq  (A7) : 


(99) 


X(s)e*‘  =  *  la  *  2-s'^Z| - 

2  _  1/2 
s  +  2s  +1 


(A15) 


Wylie  defines  the  residue  (20:818) : 


Residue  (p  )=  R(p)  =  lim  (s  -  p  )  2£(s)e 


( A16) 


The  characteristic  equation  can  be  written  in  terms  of  its  s  roots: 
s2  +  2s 1/2  +  1  =  (s1/2  +  1)  (s  1/2  +  0.5437)  (s  1/2-  0.7718  +  i  1.1151) 


*  (s  -  0.7718  -  il.  1151) 


(A17) 


The  residue  at  s  =  p  =  -0.647777+il .72127  is: 

3 

R ,  ,  =  1  im ,  ,  I .  .  . 

(pa)  «-»  (p^  L 

-1/2  <-(.  <S* 777+ il.  72127)  t  "I 

<a+C.  <5*777-11.  172127)  A  <a  +  2a  >  e _ I 

1/2  1/2  1/2  1/2  1 

<a  +  1)  (a  +  C.5437)  (a  -(.7718+11.1151)  <8  -0. 7718-H.  1151)  J 

(A18) 

The  root  at  p  can  be  written  in  terms  of  its  half-order  factors: 

9 

(S  +  0.64777-il .72127)  =  (s1/2+  0 . 77 18- i 1 . 72 127)  (s1/2+  0 . 7718+ l 1 . 72127) 


( A19) 


This  results  in  an  expression  for  the  residue: 


R.  ,  =  lim  .  , 

(P9)  (P  )3 


-1/2  .  1/2 

A  <a  +  2a  )  (a  +  O.  7718+H .  1151  )  e 


<  -O.  <S4?77+il .  72127)1 


1/2  1/2  1/2 
<8  +1)  <8  +  0.5437)  <  a  -0.7718+11.1151) 


( A20) 
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In  the  limit,  this  expression  becomes: 

R,  v  =  A  (0. 14467-iO.  1122)  e~0<s*777t  [Cos  ( 1 . 72 127t)  +  i  sin  ( 1 . 72 1 27t )  ] 

(P9} 

(A21) 

When  combined  with  the  residue  from  the  complex  conjugate  root  at 
s  =  pg=  ( -0 . 64777  + i 1 . 72 1 27)  ,  the  total  system  residue  is: 

R  =  2A  [  0.14467  cos  ( 1 . 72 127t)  +  0.1122  Sind  .721270  ]  e"0-*5*7771  (A22) 

Which  can  be  simplified  even  further  through  the  use  of  phasor  notation: 

— O .  <S4777t 

R  =  0.73232  e  cos(1.72127t  -  0.65966)  (A23) 

For  the  contour  integral  (again  refer  to  Figure  A1  for  the  paths) : 

2  2 

Along  path  2,  parameterize  s  =  <y+  io>  by  R  exp(i©)  ,  where  R  =  cr  +  go 
and  9  -  tan  *(<*>/  ©)  : 

a  -  R  exp(i©)  (A24) 

ds  =  iR  exp(i©)  d 9 

The  limits  of  integration  are  from  9  =  a,  where  a  =  tan1  (y/R)  to  n 
and  the  integral  can  be  written: 

n[R  exp(i©)  +  2R1/Zexp (-iQ/2) 3  exp (Rtexp (i©?  ]  Rexp(i©)  ^ 

O*  [R  exp(i©)]  -+  2  [R  exp(i©)]  +  1 

(A25) 

Let  the  integrand  be  denoted  3,  then: 

r i *>£  i  ,+«  >  ("I'm  jo  i 
J  I'-'  I  a-  -  j  ,v-.  , 

also 

|R  exp(i©)  |  =  |r  (cos©  +  isin©|  <  |R  |  (A27) 


(A26) 
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|R‘1/2exp(-ia/2)  J  =  |R'1X2  [cos(6M  |  ^  fl \ 


(A28) 


|exp  [R  exp(i0t) 3  | 


i 


|exp 


(Rcos0)t  ]  [  cos(Rtsin©)  + 
.  . . ^  exp  [R  cos ( t )  3 


isin  (RtsinQ)  3  |.  .  . 

( A29) 


Employing  Eqs  (A26)  -  (A29) ,  an  approximation  to  Eq  (A25)  for  large  R  is: 
A  n  D2 

”  -  J*  - - - — -  exp  [Rcos(^t)3  d&  (A30) 

zrz  a  R  +  2R  +1 

Taking  the  limit  as  R  -»  *  : 

A  r~n  M  +  op-3/2) 

lim  -  -2-  J  — ^ -  -asz  -2  exp  [Rcos(et>)  d©  (A31) 

R->  ce  2ir  a  1  t  2R  +  R 


Examine  the  exp  [Rcos(0t)3  term.  For  &  -  ny z  to  n,  cos&  <  0  and  in 
the  limit  the  term  goes  to  exp(-oo)  or  zero.  For  a  <  9  <  nszy  the  term 
is  approximately  equal  to  exp(^t).  As  R  becomes  very  large  a  £  n/2  and 
the  integral  can  be  broken  into  two  parts: 


lim 

R-»  00 


A  r  n  n/Z 

—  I  J  exp  [Rcos(0t)3  d&  +  J  exp  (H) 
L  n/2  rr/2 


(A32) 


In  the  limit,  the  integral  on  the  left  disappears  because  the  exponent 
is  negative.  The  integral  on  the  right  disappears  due  to  the  path 
length  shrinking  to  zero  in  the  limit.  Thus  the  contribution  of  path  2 
to  the  integral  is  zero.  A  similar  argument  can  be  made  for  path  6. 

On  path  4,  the  radius,  p,  is  shrunk  to  zero  in  the  limit.  Along 


1 


the  path  of  integration,  s 


cy+ico  can  be  parameterized  by  p  and  S: 


2  2 
P  =  O  + 

©  =  t  an  1  ( /  cr) 

s  =  p  exp  (  i  &) 
ds  =  ip  exp  (  i &)  d 6 


(A33) 


the  limits  of  integration  are  rc  to  (-rr)  and  the  path 


integral  is  written: 


-  A 


2rti 


p^CP  exp(iti>)  t2plxZexp(-i^/^  1  exp[  p,  exp(it»3  pexp(ifl)  id6> 
J  -  2  i/2 


[p  exp  (  i©)  ]  -+  2  C  P  exp  (  i  ^  ] 


( A.T41 


Again  let  the  integrand  be  denoted  3,  then: 


J|3|  de  >  J"|3  de| 


( A35  > 


and 

Jp  exp(i6)  |  =  | p  (cos©  +  isin£|  <  |p|  (A36) 


jp  exp(-i0/2)  |  =  |p  [cos(S/z)  +  isin(#/2)]  |  < 

|exp  tp  exp(i^t)  ]  J  =  Jexp  [  (pcos6)t  ]  t  co3(ptsin0)  + 

. . . ^  exp  [p  cos ( t) ] 

Applying  Eqs  (A35)  through  (A38) ,  Eq  (A34)  becomes: 


(P2  ♦ 


„  -1/2, 

2  P  ) 


2p 


1/2 


exp  [pcos (0t)  ]  d# 


jp"1X2|  (A37) 

isin  (  ptsin  6)  ]  |.  .  . 

( A38) 

(A39) 
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I 


In  the  limit  as  p-»  0  this  expression  becomes: 


» 

_ A 

2  TT 


-  TT 

J  0  d£?  =  0 

TT 


» 


» 


» 


and  path  4  does  not  contribute  to  the  integral. 

On  path  3,  s  is  parameterized  by  r  exp(irc),  with  the  limits 
going  from  R  to  p.  When  the  appropriate  substitutions  are  made: 


s 


=  r  exp(irt) 


=  r  (cos  n  +  isin  Tt) 


ds 

l/z 

s 

-1/2 

S 


dr  exp  (in) 

1/2 

r  exp  ( l  n/2) 
-1/2 

r  exp(*in/2) 


-dr 


ir 


1/2 


-1/2 

-ir 


r 


and  the  integral  term  is  then: 


(-r2  -  2 r  A^2) 

(  r2+  1)  +i2r1/Z 


exp 


(-rt) (-dr) 


but,  the  denominator  must  be  rationalized: 


—  J 

zm.  J 


,  2  „  “1/2.  r  ,  2  ,  ,  .  „  1^2, 

(~r  -  2r  )  [  (  r  ♦  1 )  -i2r  j 

.  .  2  ,  »  .  „  1/2,,,  2  ..  „  1/2, 
R  [(  r  +  1)  +i2r  ][(  r  +  1)  -i2r  ] 


exp  (-rt) (-dr) 


In  the  limits  as  p  ■*  0  and  R  -»  oo  : 


A 

2ra 


; 


3 

<r 


4 

r 


-1/2 

r  +  4  ♦  i2r  > 


2 

*•  2r  +  4r  +  1 


vxp  <-tV>  <h 


(A40) 


on  r 


(A41) 


(A42) 


(A43) 


(A44) 
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Dividing  through  by  i: 


A 

2TT 


3 

i.  ( r  r  +  4> 


4  2 

<r  -t-2r  +4r  +1) 


—1/2 
<2r  > 


4  2 

<r  *2r  +4r  •*•!> 


»xp  (-rt>  dr 


( A45 ) 


On  path  5: 


s 

=  r  exp(irr) 

=  r  ( 

ds 

=  dr  exp(in) 

=  -dr 

1/2 

S 

=  r  exp  (-in/*) 

=  -ir 

-1/2 

- 1  /2 

__ 

S 

=  r  exp  (-in/2) 

=  ir 

cos  n  +  isin  n)  -  -  r 

1/2 

1/2  (A46) 


and  the  limits  of  integration  are  p  to  R.  Following  the  same  procedure 
as  used  on  path  3  produces: 


A 

27T 


oof  a 

Jv.  <  r  *  r  +  4>  ^ 

4  2 

O  <r  +2r  +4r  +l> 


-1/2 
(2r  ) 


4  2 

<r  +2r  +4r  +l> 


<-rt>  dr 


( A47) 


The  limits  of  integration  on  paths  3  and  5  are  reversed.  Combining  Eqs 
(A45)  and  (A47)  and  using  limits  of  integration  from  0  to  m  produces: 


oo 


—  f  — 

n  J  4 


-1/2 
<2r  > 


®xp  <-rl>  dr 


( A48) 


O  <r  +2r  +4r  +1> 


The  total  Laplace  inversion  is  then: 


X~1i  *(s)  )  =  R  =  0.73232  e  ~°*  tf4777lco«  ( 1 .72127t  -  0.65966) 


oo 


r  _ 

n  J  4 


-1/2 

<  2r  > 


&xp  i-M)  *  (A49) 


O  <r  +2r  +  4r  + 1 > 


(105) 


The  integral  in  Eq  (A49)  does  not  have  a  simple  closed  form 
solution.  Therefore,  it  becomes  necessary  to  break  the  integral  into 


three  parts  and  make  asymptotic  approximations  to  the  upper  and  lower 

4  2 

intervals.  For  small  r,  r  +  2r  +  4r  «  1  and  Eq  (A48)  can  be  written: 


r 

2A  r  1  —1/2  „  , 

—  J  r  exp(-rt)  dr 
o 


( A50) 


The  upper  limit  of  integration  is  determined  by  a  trial  and  error 
process.  Agreement  between  the  integrands  in  Eqs  (A48)  and  (A50)  must 
be  accurate  to  within  three  significant  digits  over  the  interval  of 
integration  to  favorably  compare  with  the  experimental  data.  (Since  the 
experimental  data  contained  three  significant  figures).  The  worst  an 
asymptotic  approximation  could  do  over  the  interval  of  integration  would 
be  to  produce  a  constant  error.  The  difference  between  the  actual 
integrand  and  the  asymptotic  approximation  would  then  be  the  error  at 
the  upper  limit  of  integration  multiplied  by  the  interval  of 
integration.  If  r^is  chosen  to  be  0.00025,  then: 


-1/2 

-i/2  r 

r  =  63.426  and  - - -  =  63.1824  (A51) 

r  +2r  +  4r  +1 

and  r  =  0.00025  produces  an  estimated  error  over  the  interval  which  is 
much  smaller  than  three  significant  digits.  Eq  (A50)  contains  an 
exponential  term  which  can  be  expanded  in  an  infinite  series  and 
combined  with  the  r  1/2  term: 


-  J 

rc  J 


. 00025 

(■  .  -1/  1/2V  3/2  2 

|  (r  -r  t+r  t  +  . 

.  . )  dr  = 

2A_ 

. 00025 

r 

Z  (-1)  r <n~1/2i 

1 

o 

rr  - 

i 

o 

n  =  °  j 

( A52) 


Upon  performing  the  integration,  Eq  (A52)  becomes: 


n  n= o 


<0. 00025) 


<n+l/2>  n 


n  ♦  1/2 


(A53) 


This  series  is  convergent  for  t<4000  seconds  and  it  is  not  anticipated 
that  the  response  of  the  system  will  last  this  long. 

For  large  r,  an  integral  which  closely  approximates  the  integral 
in  Eq  (A48)  is: 


2A 

n 


oo 

r  -o/2 

J  r  exp(-rt) 


dr 


(A54) 


If  a  value  of  30  is  selected  for  r z  the  difference  between  the  integrand 
of  Eqs  ( A48)  and  (A50)  is  less  than  10  If  the  same  criteria  for 
approximating  the  lower  limit  integral  error  is  applied  to  this  upper 
integral,  then  30  is  a  good  choice.  Section  V  contains  a  discussion  of 
of  asymptotic  limit  selection.  It’s  impossible  to  perform  numerical 
integration  on  Eq  (A54) ,  but  if  the  substitution  u  =  i/r  is  made: 


1/30 

2A  r  5/'2  ,  ,  .  ,  . 

—  J  u  exp(-t/u)  du 

o 


(A55) 


and  Eq  (A50)  now  becomes: 


Using  a  composite  Simpson  integration  algorithm  (5:162-167)  the 
value  of  the  two  integral  terms  in  Eq  (A56)  can  be  determined  to  within 
a  specified  error  value.  For  example,  the  upper  integral: 


2A  r  3/2  ,  .  ,  .  , 

—  J  u  exp(-t/u)  du  = 


P  m- 1  m  "I 

-  1  f  (e|  ♦  2^  f  (x,  )  ♦  4^  f  (x„  )  +  f(i/3d  I 

3  L  s=1  ZJ=1  2J-4 


(1/30  -  O)  ,4  4 r  5/2  IL  ,  s  ,  , 

- J53—  h  flu  e*p(t/u)l|u__M  (J57) 


where 

h  s  integration  step 

/j  =  some  mean  value  (from  the  mean  value  theorem) 

2m  =  the  number  of  integration  steps  required 

{*(  )  =  the  fourth  derivative  of  the  integrand. 


The  last  term  in  Eq  (A57)  is  the  error  term.  To  obtain  a  conservative 
integration  step  for  a  wide  range  of  times,  fj  is  chosen  to  be  close  to 
the  l/aoth  end  of  the  interval  and  t  is  chosen  as  small  as  possible. 

This  selection  produces  the  largest  possible  error  for  the  interval  for 
all  time  and  appears  to  be  a  conservative  selection. 

The  same  type  of  composite  Simpson  integration  routine  can  be  used 
for  the  interval  from  0.00025  to  30.  Appendix  B  is  the  Fortran  source 
code  for  calculating  the  response  using  the  method  detailed  here.  The 


(108) 


program  is  liberally  commented  and  a  detailed  description  of  the  input 
is  included.  Using  this  program,  Figure  27  depicts  an  output  plot  of 
these  results. 


Closed  Loop  Example 

For  the  closed  loop  system  the  equation  to  Laplace  transform  is: 

x  -  1.6324  D^2(x)  +  2.891  *  -  1.1156  (x)  +  1.4142x  =  0  ( A58) 

with  inital  conditions  x(0)  =  0  and  x ( 0 )  =  A 


The  details  of  the  procedure  will  only  be  summarized  for  this  case. 
£  {  x  }  =  s*Z-  sx(  0)  -  x(0)  =  s2*  -  As 

£  {  D^2(x)}  =  s"1/2[s«2-  sx (0 }  -  x(0)  3  =  s3/2*-  As 


£  (  x  } 


=  s'#-  -  x ( 0  ) 


=  s  -  A 


( A59) 


£  {  d‘/2(x)}  =  s'1/2[s«  -  x (0)  ] 


1/2  -1/2 
=  s  2k  -  As 


Substituting  back  into  the  equation  and  solving  for  W: 


)R  = 


i/2  —1/2 

A  (3  -  1.63243  ♦  2.891  -  1.1156s  ) 

2  3/2  — 4 

s  -  1.6324s  ♦  2.891s  -  1.1156s  +  1.4142 


The  roots  of  this  characteristic  equation  in  the  s -plane  are: 


pi,pi  =  0.7892  til. 2346  ;  P2-P2=  0.0270  ±  10.8111 


Mapping  these  roots  back  into  the  s-plane  : 


Pi,Pi  =  “  °-65715  1  10.0438;  P2.P2=  -0.9014  ±  il 


.9487 


(A60) 


(A61) 


( A62) 
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The  residue  at  s  =  -0.65715  +  i0.0438  is  given  by  the  express: 


R,  ,  =  lim 
<Pt)  P 


A  (s-1 ,6324s1/2^2.891-l. 11156s"1X2)  (s±X2+0 . 027+iO ■ 81 1 1 )  exp(-st) 
(s1/2-0 . 7892  +  i  1 . 2346)  (s 1/2-0 . 7892- i 1 . 2346)  (s i/Z-0 . 027+i0 . 8 1 1 ) 


( A63 ) 


Carrying  out  the  calculations  and  combining  the  result  with  the  residue 
obtained  at  p  can  be  shown  to  produce  the  residue  expression: 

R(p  -)  =  2 . 2572 1A  exp (-0 . 65715t)  cos(0.0438t  +  0.74302)  (A64) 

The  residue  at  p  is: 

2 


.  f 

*->  p,  L 


)  exp(-st) 
iO. 811) 


CA65 ) 


A  (s~l. 6324s*X2 +2.891-1. 11156s1 ^2)  (s  t'2+0 ■ 9014- i  1 . 94877) 
(s1'^2-0.7892+il  .2346)  (s4^2-0 . 7892- i  1 . 2346)  (s  ^-0 . 027+ i 


which  results  in  the  residue: 


R,  -  =  0.82138A  exp (-0 . 90 1 4t )  cos  (1.94S7t  -  1.15255)  (A66) 


Again,  when  evaluating  the  paths  for  the  contour  integrals,  the 
only  paths  which  contribute  are  paths  3  and  5  and  the  resulting 
integral : 


_A _  j°°[2. 30854  r1/2  -  1.63425  r  1X2 1  exp  (-rt)  dr 

*  or4  -  3 . 1 1727r 3  +  7.41348r2  -  6.84149r  +  2 


(A67) 
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Finally,  after  calculation  of  the  asymptotic  approximations  the  full 
equation  in  the  time  domain: 


# 


x ( t )  =  2.25721A  exp  (  -0 . 6  57 15t )  cos(0.0438t  +  0.74302) 


+  0. 

+ 


2138A  exp ( - 0 . 90 1 4 1 )  cos 

QO  (n+l 

2.30854A  V  0.00025 
11  n=o  (  2n  +  3) 


(1 .9487t 


3) 


(-1) 


1.15255) 


+ 


_A 

n 


-0°[2. 30854  ~  1.63425  r  t''2]  exp  (-rt)  dr 

or4-  3 . 1 1727r 3  +  7.41348rZ  -  6.84149r  +  2 


+  2.30854A  r  l/2°  7^2  ....  . 

- - -  J  u  exp(-tu)  du  (A68) 

i 

When  substituted  into  the  Fortran  program,  RESPONSE,  with  the  input 
parameters  listed  in  Appendix  B,  the  output  of  this  equation  is 
illustrated  in  Figure  28. 
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Fortran  Source  Code  -  Laplace  Transform  Method 


7  MAX 
XL  I  MIT 


.  :t;°  mt-erva.  res;  :-.:e  -nc  .-.ere 

cion  ;r.terval .  From  Art  end  ix  A  API  7  ' 
AFXT'EKY/GA XLIM1T'  I '  =7:1  .  XL  I. Mil 
-2  a"  i  XL* VI”  '?  ^  -• 


:ero  a*  wr : <c r.  * o  s *  => ■ 


r*  i  n  .  «-  v  t.  _  ,  i. 

^  V  *9  9  d  1  V  „  f  *  b  V 


Lr  LAG 


eva.ua* i  or.  ox  sr.e  rest  cue  r-;: 


or..'/  wnen  3 e ' 


ooir.ts  •'  times  ■  to  evaluate 


i  r.  o  e  r  v  a . . 


AM?  ■’  I  ••  -  Amplitude  of  sinusoid  associated  with  the  ♦ 

ith  complex  conjugate  poles  * 

FREQ '  I )  -  Frequency  of  same  sinusoid  * 

PHASE (15-  Phase  shift  associated  with  the  sinusoid  * 
TAM? 'I’  -  Tamping  coefficient  associated  with  poles  * 
XMUL? f I ) -  7he  coefficient  multiplying  the  integral  * 

intervals  such  as  2 . 20254A  fir  *he  closed  * 

loop  -case  ;r.  Appendix  A  of  GA-83D-C5  * 

ITS’ I’  -  Permits  tr cublesh oo t mg  of  routines  when  * 

t.  b  ^  ^  s  ^  -i  T’  3  ^  *  O'?  *■.  b  c»  f  i  r  ^  *  t»i  rr.*  ♦ 
for  a  |:7-3r.  system.  ZiXhb  possible  outputs  ♦ 
are  available  as  follows:  * 

1 1 3 ' 1 ' :  When  set  to  0  no  debugging.  1  -pens  * 

o  h  c  o  p-  U  t,  f  *  1  * 

I D  3  i  T  i  :  Intermediate  values  for  interval  0-X1* 

IDE  '33:  Values  f  or  interval  HI  to  1.0  « 

1 3  3  >'  4  >  :  Values  for  interval  1.0  to  XT  * 

IDE!  5):  hot  used  * 

IDS  f  6 ! :  Values  for  interval  XT  to  Infinity  * 
1D3 (?) :  Not  used  * 

133(35  :  Values  from  determining  ideal  step  * 
sice  -  permits  preliminary  evaluation* 
of  program  run  time  when  N  set  to  I  * 

■  **##**»**»#*#**##**#***#«***)(*******♦♦***#***■»****♦#**♦»****«*»»»*♦*'»******* 


FILES  USED: 


INPUT . DA7  -  Input  file 
0  cTPuT • DA.  -  Contains  times  and 


■  i  a  t  e  a  r  e  s  0'  o  r.  s  e  s 
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£*»*******#***#«****♦*#*****♦**#»****<*»**»♦. «»»«»***...*♦*<•*** 

c  * 

’*  COMMON  BLOCK? 

C*  P PAP AM  -  Contains  the  information  for  evaluati r.f  res. due 

C*  response  AM?  ,  FPEQ  .  PHASE  ,  I- AX? 

C*  I  ST  I  XT  -  Ir.tegrati  on  interval  information 

c*  KL'M'T  v ~ 

•?»  H  -  Integration  step  size  for  each  interval 

C *  M  -  X umber  of  tot  al  snfcin*  ervals  .  mcludir.e  the 

C*  er. dp:  i nts  ,  m  each  interval 

C*  I 3  CO  XT  -  Debugging  control  parameters  ILB'l-? 

C  * 

C#******«**##***«****#*#«****<**<*****»****f<*******  -»********♦*♦**»»■***♦*>!**»' 

C  use  SDEBUG  t:  help  identify  reasons  for  program  failure 
SDE3UC- 

:m?l icit  locbli  precision  ja-h.o-z) 

COMMON  /  P.F  AKA  M/  A  Mu'  t  2 )  ,  FP.EQ  ( 2 !  ,  PHASE  !2)  ,  LAMP  ( 2  ! 

COMMON  '  1  NT  I  NT -■  XLIMIT'3!  , XMULT ! 3 ;  ,H!3’  .PI , M ( 3  5 
COMMON  ■ D3CCNT  '  I  LB 'S' 

?:=3.  141" 9 1 554 

("*  ,  -  •**  ^  *  »•».  *•*,«<  •*  a  **  y*  -  .,«•*■...*•  ^  ^  ■? 

m*r£  ,  i  -  ,  i  V  tilU  V  a*  b  •* 

OPEN  • 5 . FI LE - ' IN? L'T  .  OAT '  ,  STATES  =  ‘  OIL '  . 

OPEN  w.FIlE^'OVTPTT.LAT'  ,3?A7*JS= ’NEW” 

/*» 

C  read  inrut  raraxeters 


PEAn 

re  *  \  •"C'f  At" 

■  v  (  "  -  .  rt.*  .  , 

TTEPO , 7MAX 

HEAP 

i  5  ,  * ;  c AM? ' I ? 

, FREQ ( I )  ,  ? 

PEAL 

•?,»)  '  XMULT  l 

: >  .i=!,3) 

PEAL 

' 5 , ‘ )  N . ( IL3 ( 

n  .3) 

intermediate  values  are  to  be  printed  open  debug  out pm 
IF  fIDBC)  .  EQ.  1)  OPEN  (3  ,  FILE  -  ’  DEBUG.  DAT '  ,  STATUS  =  ’  NEW 
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1*1  UJ 


(121) 
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call  the  function  routines  to  evaluate  A.  2 ,  and  -  and  the  lire: 
fourth  derivatives  A1 .  A2 ,  A3.  A4  ,  ev: 

A* AC Aid  K2AS,Ni 
A1=A?CALC  .'ZBAS.X: 

A2  =  ACC ALC  ;  NEAR ,  X  ’ 

A?  =  AT  CMC  ■  X3AB .  X ) 

A  4  =  A  v  C  A  L  C  '  X3  AP. ,  X  •• 

2  =  HCALC i MPAS ,  7 ,  X! 

B1=H?CALC;X3AX.T,X) 

32;  =  33CALC (XBAS , T , Xi 
33 -BTC ALC ' XBAS , T , X) 

34=3QCALC iXEAH.T.Xi 
0=3 CALC • XBAR.X! 

C  1  =C?CALC ( XBAS , X 1 
C 2 -CSC ALC : XBAS . X! 

T  3  =  CTC ALC  f  XBAS . X ' 

C4  =CQCALC f XBAS , X ’ 

the  chain  rule  is  used  t:  combine  the  A,  3,  ar.d  C  derivatives  in; 
fourth  derivative  for  the  expression  •' AB/C’ 


72=  f  1  .  27i>A2*B»  *C 
- ! 1 . 2 31 »A2 *3 


*  ♦  2  5  /  .  C « *  3  ;  i-  5 . i>I'0*A2*3*C2,f  |C**2'  i 
r  c.  1  /  :  C » ♦  2  ■  r  -  .  6  .  C  C  C  ♦  A  2  *  B  2  /  C  ’ 


-  '  -  .  0D0*  A1  *3*C3  '  C*  *2 ' 


73=  •;  2  .  431*A1*3»C1*C2  '  ;C**3  :  •  -  ;  2 . 431*A1 *2»  'Cl*«2:  ■  C»* 

-  '4 . 0D0»A1*3*C3-'"C**2'  1  *  '  2. 431 *A1 *3 1  *  1 C 1  *  *  2 ) / ' C  * 

;  -  ;  1  .  2  3 1  *  A 1  *  B 1  *  C  2  /  <  C  *  *  2  ;  ;  -  '  1  .  1 D 1  *  A 1  *  B  2  ♦  C  1  /  ( C  *  *  1 i  ' 

?  *  '  4  *  A 1  *  3  3  /  C ! 

74=  ;  2 . 4D 1  *  A*  3*  ( C 1  *  *4 )  -  <C«  *5 )  1  -  f  3  t  L 1  *  A*  3*  (C 1  *  *  2 )  *C2  •'  < 
*  ;  &  .  CD0*A*3*C  1  *C3/  (0*3)  '  -  !6. ODO*A*3*  (C2**2  •  /  iC* 
:  -  :  A*3  »C4 /  ! C* »2  >  !  *  i A *34  C )  -  1 4 . OLC  * A*33  »C  1  /  u:>»2  : 

7 *  =  >'  2  .  4  T  2  *  A  *  3 1  *  C 1  *  C  2  •' !  C  *  *  3  )  1  -  f  2 . 4  7 1  *  A  *  3 1  *  C 1  *  *  2  '  •'  :  C  *  * 

-  :  4  .  c:.0«A*31  *02/  <  ) 

76= ; A 4  *3  C i - ; 6 . 030* A*B2*C2/  C  +  *2'  : * •  1  . 2D1 » A  *32 *  f  C 1  *  *  1 
?QTAD  =  7 1*72*73*74+75+76 


2**4)  ) 
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O.  Vi. 


(126) 


05  hi 


fjJ  »ii  l»j  hj  *  *  *  i  it>  »:*.  **  «'  »  i  O  r*J  p.) 


r>5  i*j 


DO  VE1E  7  HE? I  S'  I  ?S  EVICTION  CCAEC  f  XPA7; .  X ' 

: m?l :c :t  :-ovbl:  ?3ec:s:ox  (A-h,o-z> 


•his  function  ev-aLuases  der.onir.anor  :c ' 


tec  interned  lane  va:u 


terva.  .  r.  o  iru.  oan-c 


:  r.  ter  v  a .  ■  _i .  t 


indicates  tr.e  ucper  .  nterva.  '  ~:2  tc  infinity  i s  to  be  eva.ua  ted 
VOTE:  when  changing  tc  a  different  1  -..near  esuaticn  this  :u.tct_ 

n--s  t  be  charged  to  accurately  reflect  the  integrand  :  e : .-. 
evaluated 

IE  ,  L T .  3  ;  THE',' 

OCAI.C  -  XEA3**  4.c*2. 0  TO  »  EE. A?.*  *2  .  c*4 . 01 C  ♦MEAT.*  *  1 .  c*X3A?:*  *0 

tt  c; 

COAL C=  I. EDO 


TO  TELE  ?  EE  CIO  ICE  EVNCTI  IE  OEOALO  'EE  AH  .  M  ■ 

I  VET.  I  T  IT  M  TEI-E  E-EEC-  ITT  OX  '  A  -H  ,  1  -Z 

this  function  evaluates  the  first  derivative  o:  the  denominator  plyr.rr.: ; 
cf  the  '.  r.  t  e,=  r  a  r.  d  at  a  selected  ir.  termed i ate  value  ir.  the  interval, 
n  3  indicates  the  middle  interval  El  tc  EE'  integrand  is  t:  be  evaluate: 
r  =  3  indicates  the  upper  interval  -El.  tc  infinity  is  tc  be  evaluated: 

VOTE:  when  changing  to  a  different  linear  equation  this  fun: tier. 

mt:s t  he  changed  to  accurately  reflect  the  .ntegrand  being 
evaluated 

IF  fX.LT.?’  THEN 

Or TALC=  4 . 5D0«X3AS**3 . 5*5 . 0D0«X5AS**  1  .  t-c  .  0T0*M2A?.«*0  .  t 
*5  .  OD- 1 *XEAE* * '  -  0 .  t : 


C?CALC=  O.CDO 


stv*n  t  r 
n  ^  i.  r 


RETURN 


(131) 


(134) 
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Appendix  C 

Mi ttag  2.  Lef  f ler  Expansion  Method  f or  Predicting  Time  Response 


The  Laplace  inversion  method  for  analytical  prediction  of  closed- 
and  open-loop  system  response  agrees  closely  with  experimental  results. 
The  method  has  drawbacks  which  ma  e  alternative  methods  requiring  less 
up-front  work  attractive.  Up-front  work  refers  to  hand  calculation  of 
the  residues,  calculation  of  four  derivatives  of  three  functions,  then 
programming  these  derivatives  into  the  function  routines,  and  finally 
making  engineering  judgments  on  the  points  at  which  to  apply  the 
asymptotic  limit  expansions.  There  are  many  points  at  which  errors  can 
be  made  in  this  process.  Is  there  a  simpler  analytic  tool  to  use’ 

The  Mittag-Lef f ler  expansion  method  (12)  employs  a  time-domain 

3/2  1  1/2 

solution  technique.  Identity  relations  for  (x) ,  D^tx),  and  Dt  (x) 
exist  as  follows: 


and 


,  ..i/2 

Dt  (x)  =  Dt  x 

Dt  [D t  (x)  ]  =  D^x  = 

-t/2,  1,  .  ,  ,.3/2 

Dt  [Dt  (x)  ]  =  Dt  x 

,.1/2,  3/2,  ,  .  _  2 

Dt  [D  t  (x)  1  =  Dtx  = 


(Cl) 


The  entire  second-order  equation  can  then  be  written: 

mD*/2[D^/2(x)  ]  +  ctD  J/2[D  J(x)  ]  +  c  £  *^x)  ]  +  c  D  *^fx)  =  kx  =  u(t) 

(C2) 
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Eqs  (Cl)  and  (C2)  can  be  put  together  into  a  set  of  linear  equations  in 
matrix  form: 


„l/2  ,  % 

Dt  (x) 


O  O  O  1 

O  O  1  o 

O  1  o  o 

m  c  c  c 


r  rs3y2  > 

Dt  x 

Dt  x 
- 1/2 

Dt  x 

-  + 

.  x  > 

OO-l  o 
0-100 
-1  o  o  o 
O  O  O  k 


.1/2  , 


f  „  3/2 

D  t  x 

c 

D  t  x 

o 

•i 

Dt/Zx 

-  =  * 

c 

x  . 

u<t ) 

k.  / 

f  C  3 ) 


If  the  fractional  derivative  operator,  (x)  ,  is  replaced  by  X, 
Eq  C3  becomes: 


o 

o 

o 

m 


o 

o 

1 


O  1 

1  o 
o  o 


C  c  c  . 

1  2  3-* 


0  0-1  O 

0-100 
-1  o  o  o 

O  O  O  k 


> 

fx3  1 
x2 

o 

o 

‘  * 

X 

-  ~  H 

o 

-1  j 

1 

V.  > 

u<t> 

t  / 

(C4 ) 


and  if  u(t)  =  0,  the  zero  state  response  of  the  system  is 
characterized : 


o 

O 

o 

-1 

-l 

X 

1 - 

^  0 

f  ^ 

X2 

f 

\ 

o 

o 

-  1 

X 

o 

o 

1 

X 

•  —  ■* 

o 

mX 

c  X 

l 

c  X 

z 

c  X+k 
3  J 

1 1 J 

o 

(C5) 


Now  this  is  an  eigenvalue  problem  and  has  the  same  characteristic 


1/2 


equation  of  the  system  posed  in  terms  of  Dt  : 


4  3  2 

mX  +  c  X  +  c  X  +  c  X  +  k 

1  2  3 


(C6) 


Using  a  root  solving  routine,  such  as  that  contained  in  Matrix^  ,  the 
eigenvalues  are  simply  the  roots  of  the  characteristic  open-  or  closed- 
loop  differential  equation.  The  system  response  can  now  be  posed  in 
eigenvalue/eigenvector  format: 
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3/2  " 

D  x  ( t) 
t 

Xi 

xi 

xi 

x!  ' 

r  *v 

a  (t) 

i 

D1  x(t) 

t 

-  = 

Xi 

xl 

xi 

xi 

a2(t) 

D^x(t) 

Xi 

X2 

\3 

X A 

a3(t) 

.  x(t)  . 

1 

1 

1 

1 

.a4(t)  . 

or  symbolically: 

{  D (t)  }  =  [  A  ]  {  A(t)  ) 

Total  system  response  is  the  sum  of  the  characteristic  behavior 
belonging  to  each  root  of  the  differential  equation: 

4 

x(t)  -  Z  a.(t) 

t  =  1  V 

and  the  a  (t)  is  a  series  expansion: 

a  (t)  =  a  (0)  E  (  Xt 

X.  L  1X2  V. 


The  expansion,  is  defined: 


oo 

E  (z)  =  £ 

1/2  n=0 


ru  +  n/2) 


r  is  the  gamma  function  or  generalized  factorial.  The  total 
solution  is  then: 


The  aJO) 
Eq  (C8)  : 


x  (t) 

terms  arise 


,  .  1/2,  n 

(  Vt  ) 


=  Z 

L  —  O 


Z 

n  =  0 


a.  (0) 


r 


from  the  initial 


(  1  +  n/2  ) 
condition  vector, 


{  A ( 0 )  }  =  [  A  ]~l  (  D (0)  } 


{  D (0)  }. 


The  a^(t)  terms  can  be  complex  numbers.  The  solution  can  be 
implemented  on  a  mainframe  computer  using  any  type  of  programming 
language.  However,  to  implement  it  on  a  Z248  PC,  one  must  use  the 


(C7 ) 

(C8) 

(C9) 

(CIO) 

(Cl  1) 

(C12) 

From 

(C 13) 
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macro-programming  routines  in  Matrix^  to  be  able  to  work  with  complex 
numbers.  These  macro-routines  are  portable  to  the  Matrix^  version 
contained  on  the  AFIT  Vax  mainframes.  Appendix  D  contains  the 
Matrix^  routines,  files,  and  instructions  for  their  use.  It  is 
programmed  so  the  expansion  E  can  be  done  for  any  i/nth  order  value 
desired.  This  is  accomplished  by  responding  to  the  prompt  for  system 
dimension  with  (2*n) .  For  example,  if  the  second  order  equation 
involved  quarter  derivatives  instead  of  halves  the  response  would  be 
[  (2)  (4)  ]  =  8. 


(139) 


ft 


$ 

Matrix 
- x 


Appendix  D 

Code  for  Mittag  ^  Lef f ler  Expansion 


Page  141  contains  the  code  for  MITLEFR.DAT  which  is  an  executable  file. 
To  run,  execute  Matrix^,  ensure  MITLEFA.DAT  (page  142)  and  MITLEFB.DAT 
(page  143)  are  present  on  disk,  then  enter  ‘EXECUTE  [ 'MITLEFR.DAT' ] ‘ . 
Follow  the  prompts. 


(140) 


(142) 


// [TOUT ,XOUT , XERR1=HSPCAL ( AZEEO , A,I) 

II... 

II ...  accept  routine  control  input 

II... 

INDI1S  TSTABT  ’IIPDT  STABT  THE:  ' 

INDI1E  TZEBO  ’IIPDT  STABT  TIME  OF  DESIRED  IITEBVAL:  ' 

INQUIRE  TMAX  ’ INPUT  MAX  TIME:  ’ 

INQUIRE  MUMPTS  'IIPDT  SUM  OF  POINTS  TO  CALCULATE:  ’ 

INQUIRE  IUMITS  'MOM  OF  TEBMS  II  MITLEF  SERIES:  ’ 

II... 

II...  calculate  tine  increment  A  setup  to  correctly  define  initial  time 

II... 

TIN*  (TMAX-  TZEBO) /IDMPTS; 

T=TZEBO-TIMC; 

FOB  IK«0: IDMPTS,... 

II ... 

II  establish  matrixx  compatible  storage  index  and  initial  increment  values 

//... 

IIDX-IIC+1; . . . 

ZX=0; . . . 

W+TIIC;... 

II ... 

II  if  Interval  start  prior  to  response  start  set  oalcuaition  time  to  zero. 

II... 

IF  T<TSTABT ,TCALC=0 ; ELSE  TCALC=T-TSTART;EHD , . . . 

II... 

II  set  up  initial  values  for  gamma  function... 

II... 

PRE1= 1 . 0; PRS2=0 . 5 ;LBIT= 1 ; . . . 

II ... 

II  calculate  the  first  2  terms  in  the  Mittag-Leffler  (ML)  series... 

//  ...  use  these  terms  as  multipliers  of  the  integer  and... 

H  non-integer  gamma  function  terms  of  the  M-L  series... 

II... 

FOB  1=1:1,... 

Z(I) =-A(I) *SQRT(TCALC) ; . . . 

ELASTHIM.0/PBE1;... 

ELAST2 ( I) =Z ( I) / (SQBT (PI ) *PRE2) ; . . . 

ZX= ZX+ AZERO ( I ) * ( ELAST 1 ( I ) ♦ ELAST2 ( I ) ) ; END _ 

// . . . 

//  set  up  the  correct  denominator  for  the  next  non-integer  gamma  tern... 

II ... 

PRE2=PRE2+ 1 
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II... 

//  L8IT  controls  the  use  of  the  integer  or  non-integer  gamma  function  terms... 
//  ...the  next  loop  calculates  the  j+2  term  of  the  M-L  series  for  each  of... 

//  the  n  roots. . . 

//. . . 

FOB  J=1:MUMITS,... 

B=0+JAY*0;... 

II ... 

//  calculate  the  gamma  function  term  and  set  up  for  the  next  loop... 

II ... 

IF  LBIT=1 ,DEM=PRE1 ;PHE1=PRE1+1 ;LBIT=2; . . . 

ELSE  DEI=  PHE2 ; PEE2= P8E2 ♦ 1 ; LBIT= 1 ; . . . 

END,... 

FOB  1=1:1,... 

II ... 

II  multiply  the  last  integer/non-integer  term  by  z»*2/n-l  to  get  the... 

//  correct  M-L  term  for  the  ith  root... 

II ... 

IF  LBIT=2.  ELAST1 (I)=(ELAST1 (I)*2(I)»*2) /DEM; . . . 
B=B+A2EB0(I)*ELAST1(I):... 

ELSE  EL AST 2 ( I )  =  ( ELAST2 ( I ) *Z ( I ) « *  2 ) / DEM ; . . . 
R=BUZEB0(I)«ELAST2m:... 

END,... 

EHD .... 

II ... 

II  ZX  represents  the  sum  of  the  l  roots  for  each  term. . . 

II... 

ZX=ZX+B; . . . 

ESD,... 

//. . . 

II  store  the  output  for  each  time  step... 

t  I 

XOOT  (IMDX)=BEAL(ZX) ;XEBB(IMDX)=IMAG(ZX) ;T0OT(IIDX) =T; . . . 

EHD 

II... 

II  save  the  data  file  in  a  format  readable  by  grapher... 

II... 

FSAVE  ’MITAPP.DAT’  TOOT  XOOT 
RETF 
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A  summary  of  the  design  and  fabrication  process  involved  in 
realizing  the  Oldham-Zoski  resistor/capacitor  domino- ladder  circuit, 

V  x> 

which  will  perform  the  operation  d/dt  on  an  input  electrical  signal 
is  presented.  The  valid  range  for  v. 

-1<  v  <  l  (El) 

where  v  less  than  zero  indicates  integration. 

The  basic  circuit  is  referred  to  as  a  dorm.no  ladder  and 
consists  of  a  chain  of  resistors  and  capacitors,  the  two 
chains  being  connected  at  each  node  as  in  Figure  3.  Each 
resistor  is  a  constant  factor  multiple  of  its  predecessor 
as  is  each  capacitor  giving  the  relationship: 

R  S  R  g”J  and  c.  =  C  a  S  (E2) 

J  O  *  JO 

where  both  g  and  a  are  greater  than  unity. 

To  proceed  select  a  value  of  V>;  then 

In  a  <  3/2  t>2/3  (E3) 

and  In  g  =  <1— V»>  V  *  In  a  <l?:3l>  (E4) 

The  rest  of  the  method  follows  the  basic  guidelines  established  in 
(17).  However,  the  derivation  is  put  in  terms  of  frequency  instead  of 
time,  as  was  done  in  the  original  Oldham-Zoski  article.  First,  a  minimum 
frequency,  f  is  selected  for  the  application  (fmis  in  Hz) .  The  time 
constant  of  the  first  resistor-capacitor  pair  is  then: 


H  C  =  HI  exp  (,.-.^2X3) 

o  o 


sec 


(E5) 


*  Gg 

m 

Any  combination  of  resistors  and  capacitors  which  produces  this  time 
constant  is  acceptable.  As  explained  in  the  analysis  section,  the 
smaller  the  capacitance  values  the  better.  This  must  be  balanced,  on  the 
other  hand,  by  the  values  of  the  largest  resistors  and  smallest 
capacitors  available.  If  |l*|  is  =  1  <  or  the  frequency  range  is  large, 
the  design  will  require  more  cells  than  if  |i>|  is  small.  This  could 
possibly  put  a  limit  on  the  initial  resistance  and  capacitance  values 
used.  It  is  also  wise  to  make  the  tolerance  on  the  components  less 
than  IX.  This  appears  to  guarantee  good  performance. 

To  calculate  the  number  of  cells  required  in  the  domino-ladder, 
select  a  desired  upper  frequency  limit,  f^.  The  number  of  cells 


required,  N,  is  then: 


2/9  —  1 

N  +  1  >  [  5.5  +  ln(  i'i  )  -  3w  ]  [ In (gG)  ] 

M  m 


(E6) 


Some  enhancements  for  the  high-  and  low-frequency  performance  are 
available  (17:33).  To  increase  the  accuracy  at  high  frequencies,  modify 
the  final  resistor  capacitor  pair: 


(E7) 


• 

from 

Bn  *  Ro*‘N 

to 

1X2  R^-N 

c  =  c  g"n 

N  O 

to 

2  C  G_N 

O 

• 

and  introducing  a  final 

resistor  at 

the 

output  of  t 

Bm.,  -  Bo  «' 

N(ln 

g)_1 

(E8) 
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Low-frequency  performance  can  be  improved  by  adding  an  additional 
resistor-capacitor  pair  at  the  start  of  the  circuit  : 

{  C  G 
R  =  (  1/2  *  — -)  R  t  and  C - - -  (E9) 

-1  In  <3  Cr  -1 

[  1/2  +  1/  InO  ] 

While  circuits  produced  using  this  method  possess  the  correct  Bode 
magnitude  slope  and  phase  response,  the  actual  gain  value  at  u  =  1 
requires  adjustment  to  0  dB.  This  adjustment  depends  on: 

v  s  fractional  order 

C  =  Capacitance  of  the  analog  circuit  interface  component 

R  =  Resistance  of  the  analog  circuit  interface  component 

R  =  Resistance  of  the  Oldham-Zoski  base  resistor 

o 

Cq  s  Capacitance  of  the  Oldham-Zoski  base  capacitor 

On  the  Pace  analog  computer  (model  TR48)  there  are  two  possible 
interface  resistor  values,  10KO  and  100KO.  The  two  interface  capacitor 
values  are  0.02AiF  and  lO.O^F.  Figure  47  shows  the  possible 
configurations  available  to  obtain  specific  fractional  differentiator  or 
integrator  values  and  the  corresponding  equations  to  calculate  the  gains 
required  for  a  given  circuit.  Table  2  contains  the  calculations  of  the 
gains  for  the  two  circuits  built  for  this  thesis  and  compares  predicted 
and  actual  gains  required.  Table  3  details  the  method  of  configuring 
the  operational  amplifiers  and  capacitors  to  take  advantage  of  the  four 
different  interface  components . For  further  information  on  the 
Oldham-Zoski  circuit  refer  to  reference  17. 
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Table  2.  Circuit  Gains  -  Predicted  Versus  Experimental 
Oldham- Zoski --Analog  Interface  Con f i jurat i on  4 


IBSSHi 

■£9 

Predicted 

iti  t: 

Experimental 

Circuit  1 

28 . 52 

25.08 

Circuit  2 

13.28 

11.47 

Table  3.  Details  for  Realizing  the  Four  Different  Analog  Computer 
Interface  Configurations  for  the  Qldham-Zoski  Circuit 


CONFIGURATION  1 


Input  signal  to  either  a  xlO  terminal  on  an 
operational  amplifier  (R=10k)  or  a  xl  input 
terminal  (R=100k).  Output  of  the  amplifier 
is  input  to  the  OZ  circuit,  output  of  the 
OZ  circuit  is  connected  to  the  xlO  input 
terminal  on  the  amplifier  (10X0  must  removed 
from  the  last  OZ  resistor  before  connecting. 


CONFIGURATION  2 


CONFIGURATION  3 


CONFIGURATION  4 


Remove  T  plug  from  one  of  the  integrator 
modules.  Insert  Z  pin  bottle  plugs  in  the 
* 0 .  1(3'  and  'SPEC*  terminals,  and  a  4  pin 
bottle  plug  in  the  'OPR/RESET'  bus.  In 
this  configuration  C*10^.  Removing  the 
bottle  plug  from  the  ‘0.10 (T  area  gives 
C*lwF.  The  OZ  circuit  is  connected  to  the 
amplifier  as  in  Configuration  1.  The  input 
to  the  capacitor  is  connected  to  the  'SJ' 
input  terminal,  the  output  is  then  avail¬ 
able  at  the  'O'  terminal.  Connect  this 
terminal  to  an  input  terminal  on  the  oper¬ 
ational  amplifier.  (Jse  only  a  2  pin  bottle 
plug  on  the  operational  amplifier  to  con¬ 
nect  the  'B'  and  'SJ'  terminals. 


Configure  the  integrator  as  in  Configura¬ 
tion  2  with  the  following  exception  -  the 
input  to  the  integrator  (capacitor)  is  the 
output  signal  of  the  operational  amplifier 
The  input  signal  is  connected  to  the  input 
terminal  of  the  OZ  circuit.  The  output 
terminal  of  the  OZ  circuit  is  then  con¬ 
nected  to  the  input  of  the  amplifier. 


Configure  the  operational  amplifier  with 
a  4  pin  bottle  plug  (RslOOKii.  If  R=10K 
is  desired,  use  two  2  pin  bottle  plugs  - 
one  connecting  "SJ*  and  'B'  and  one  con¬ 
necting  an  output  terminal  to  a  xlO  input 
terminal.  The  input  signal  is  connected 
to  the  input  terminal  of  the  OZ  circuit. 
The  output  terminal  of  the  OZ  circuit  is 
connected  to  the  input  of  the  amplifier. 
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Detailed  Experimental  Procedure 


For  each  task  listed  in  Section  IV  this  appendix  details  the 
instrumentation  configuration  used  to  experimentally  gather  the  data. 

It  also  provides  references  to  equipment  manuals  should  additional 
information  concerning  the  equipment  be  of  interest. 

Circuit  Build 

The  Tektronix  Digital  Multimeter  (model  DM501) (28)  and  Dynascan 
Capacitance  Meter  (model  820)  are  required  for  this  step.  The  scale 
on  the  multimeter  should  be  set  to  ohms  and  the  proper  range  selected 
for  the  resistor  being  measured.  A  note  of  caution  on  the  capacitance 
meter:  as  the  capacitances  get  smaller,  the  offset  error  becomes  large. 
Take  this  into  consideration  when  measuring  pF  capacitors. 

Circui t  Test 

Equipment  Used 

The  criteria  for  a  fractional  derivative  of  order  one-half  is 
a  Bode  plot  with  a  slope  of  lOdB/decade  and  a  constant  phase  shift  of  45 
degrees.  To  determine  if  this  criteria  was  satisfied  by  the  circuits 
fabricated,  a  function  generator,  counter/timer ,  and  oscilloscope  were 


Tektronix  Function  Generator  (model  FG506)  (24:1-1,1-3). 

Capable  of  producing  a  10  volt  peak-to-peak  sine  wave  at  frequencies  from 
.001  to  50  MHz  when  shunted  across  a  500  impedance  load. 

Tektronix  Universal  Counter/Timer  (model  DC509)  (25:2-1,2-13). 
Two  channel  input,  with  adjustable  trigger  voltage  for  each  channel. 
Capable  of  measuring  frequency  to  within  1x10  a  Hz  for  the  frequencies 

-9 

of  interest;  period  to  within  1X10  sec;  and  phase  difference  between 

-JS 

channels  to  within  1x10  sec  when  triggers  are  set  correctly. 

Tektronix  Oscilloscope  (model  SC504) ( 26 : 1 . 1 -2 . 29) .  Two  channel 
display,  graticule  on  the  screen  capable  of  representing  0-100%  of  5 
times  the  magnitude  of  the  volts/division  setting.  Resolution  is  ±2%  of 
this  range.  Volts/division  settings  spanning  5mV  to  10V.  The  unit  has 
adjustable  horizontal  sweep  frequencies. 

Using  BNC  T  fittings  and  coax  cables,  the  function  generator 
output  is  connected  to  Channel  A  of  the  counter/timer ,  Channel  1  of  the 
oscilloscope,  and  the  input  of  the  Oldham-Zoski  circuit.  This  circuit 
produces  an  output  current  which  is  the  fractional  derivative  of 
the  input  voltage.  The  voltage  across  the  last  resistor  is  just  the 
value  of  this  current  multiplied  by  the  value  of  the  resistor.  This 
voltage  is  the  voltage  to  be  compared  with  the  input  voltage.  (It  is 
important  to  place  the  output  measurement  device  in  parallel  with  this 
last  resistor.  Since  they  are  high-impedance  devices,  very  little 
current  will  be  required  to  produce  the  measurement,  whereas  a  series 
measurement  would  seriously  affect  the  circuit  current).  The  voltage 


across  the  last  resistor  is  connected  to  Channel  B  of  the  counter/timer 


and  Channel  2  of  the  oscilloscope. 


Procedure 


The  gain  of  the  circuit  is  defined,  in  dB,  as: 

20  log  (v  /  v  )  (FI) 

lO  out  m 

It  is  not  necessary  to  maintain  the  same  input  amplitude  over  the  range 
of  frequencies  tested.  However,  it  leads  to  the  conclusion  that  the 
voltage  input  to  the  circuit  is  not  frequency  dependent.  Since  the 
performance  of  the  circuit  is  necessarily  dependent  on  frequency,  it 
makes  the  calculation  of  gain  straight  forward  if  the  input  voltage  is 


constant . 


Without  the  circuit  attached  to  the  function  generator,  a  500 
resistor  should  be  shunted  across  the  output  of  the  function  generator, 
and  then  a  common  connection  should  be  made  to  ground;  this  matches  the 
impedance  of  the  function  generator.  A  sweep  across  the  frequency  band 
to  be  investigated  (0.01  to  1000  Hz)should  be  made.  The  amplitude  of 
the  input  voltage  on  Channel  1  of  the  oscilloscope  should  be  observed. 

If  no  changes  are  evident  then  the  assumption  can  be  made  that  the  0Z 
circuit  input  voltage  remains  constant.  The  circuit  should  be  connected 
in  parallel  with  the  500  resistor.  Since  the  impedance  of  the 
Oldham-Zoski  circuit  is  several  orders  of  magnitude  higher  than  the  500 
resistor,  the  parallel  combination  will  appear  as  a  500  impedance  to 


the  function  generator. 


Amplitude  Ratio  Determination 

The  procedure  is  as  follows: 

1.  Set  function  generator  to  a  desired  test  frequency  (record) 

2.  Set  counter/timer  to  frequency  and  validate  the  frequency. 

3.  Ground  both  channels  of  the  oscilloscope  and  adjust  the  zero 
position  of  each  channel. 

4.  Set  both  channels  to  DC  input  coupling 

5.  Adjust  the  volts/division  reference  for  each  channel  to  give  the 
maximum  deflection  without  exceeding  the  1 00'/.  mark  on  the 
graticule . 

6.  If  a  volts/division  change  was  made,  repeat  step  2  for  the 
appropriate  channel. 

7.  Use  the  chop  display  mode  to  look  at  both  signals  at  once. 

8.  Use  the  percent  values  listed  on  the  left  side  of  the  graticule 
to  record  each  channel's  amplitude  as  a  given  percentage  of  5x 
the  value  of  the  volts/division  setting  for  that  channel. 
Divisions  on  the  left  hand  side  of  the  graticule  are  4%  per 
minor  division,  giving  accuracies  of  ±27.  Use  Eq  (FI) 

to  determine  the  gain  for  the  frequency  just  tested. 


Phase  Determination 


0 


0 


0 


After  calculating  the  gain  at  a  given  frequency  the  following 
procedure  will  determine  the  phase  angle: 

1.  On  the  counter/timer  select  the  ‘PERIOD  A’  function. 

2.  Select  the  desired  number  of  cycles  to  be  averaged  together  on 
the  ‘Averages'  selector.  Select  1  for  frequencies  less  than 
0.1  Hz,  ‘Auto'  for  frequencies  up  to  1  Hz,  and  adjust 

to  stabilize  the  reading  for  higher  frequencies.  Record  the 
period  in  the  data  table  being  made. 

3.  Adjust  the  trigger  voltages  to  provide  an  accurate  evaluator 
of  the  phase  shift. 

3a)  For  both  channels,  ensure  the  following  configuration: 
‘SLOPE*:  both  the  same 
‘ATTEN':  *xl* 

‘SOURCE':  'Ext" 

‘ COUPL' :  AC 

Ensure  the  button  on  the  timer/counter  marked  ‘AUTO  TRIG 
LEVEL"  is  not  depressed. 

3b)  Using  the  multimeter  with  needle  probes,  select 

*DC  VOLTS'  and  ‘2 00  mV',  insert  the  positive  lead  into 
the  contact  marked  ‘A  -  TRIG  LEVEL’  on  the  timer/counter 
and  the  negative  lead  into  the  contact  marked  ‘COM  - 
SHAPED  OUT*.  Adjust  the  channel  A  trigger  level  to  as 
close  to  zero  as  possible  on  the  multimeter  readout. 
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Switch  the  positive  lead  to  "B  -  TRIG  LEVEL"  and  repeat 
the  adjustment  for  channel  B  to  get  the  trigger  as  close 
as  possible  to  the  value  for  channel  A. 

4)  Select  the  "TIME  A-*B “  function,  the  phase  shift  time  will 
appear  in  the  readout.  Record. 

5)  Calculate  the  phase  shift  as  follows: 

("TIME  A-»B'/  Period)  *  360  =  <t>  (degrees)  (F2) 

6)  Record  values;  repeat  gain  and  phase  procedure  as  required. 

Hewlett-Packard  HP7090A  Analog  Measurement  System 

Overall  circuit  performance  depends  not  only  on  steady-state 
performance,  but  also  on  start-up  transient  characteristics.  It  became 
clear  that  observing  the  startup  transient  was  a  necessity.  How  long 
did  it  take  the  circuit  to  become  an  efficient  fractional 
differentiator’ 

To  answer  this  question,  a  Hewlett-Packard  Analog-Digital 
Measurement  Recording  System  (model  7090A)  was  obtained  and  substituted 
for  the  oscilloscope.  The  recorder  is  capable  of  plotting  3  channels  of 
input  data,  with  a  sampling  rate  of  33  kHz.  The  model  7090A  also  has 
three  1000  word-length  buffers  (one  buffer  for  each  channel).  Each  word 
represents  a  time  slice  of  1000th  of  the  total  time  window  specified 
prior  to  data  recording.  After  experimenting  with  this  total  time 
window,  it  was  observed  that  three  to  five  cycles  of  the  input  sine  wave 


0  produced  a  steady-state  response  of  the  circuit,  and  sufficient  data  to 

accurately  determine  gain  and  phase  shift.  The  equation  used  as  a  guide 
for  the  total  time  setting: 

•  t  =  5  cycles  *  C  1/  (  f  ) ]  (F3) 

total.  Hz 

Another  feature  of  the  model  7090A  is  the  ability  to  record  data  prior 
to  the  start-up  of  the  system.  A  trigger  function  represents  the  actual 

0  start  of  the  measurement.  This  trigger  can  be  set  internally  to  occur 

at  a  specific  clock  time,  or  it  can  be  controlled  externally.  By 
specifying  a  duration  of  pre-trigger  recording  time,  t  ,  all  events 

pt 

0  occurring  within  t  seconds  prior  to  the  actual  trigger  event  are 

recorded.  In  addition,  if  t  seconds  pass  after  starting  to  fill  the 
buffers  without  the  trigger  event  occurring,  the  oldest  pre-trigger  data 

•  is  written  over  so  that  only  the  latest  t  seconds  of  data  are  retained 

pt 

in  the  buffer..  These  features  allow  the  actual  startup  transient  of 

the  circuit  under  test  to  be  recorded  (27:3-1,24). 

I 

0  The  external  triggering  system  of  the  model  7090A  can  be  accessed 

through  a  BNC  fitting  on  the  right-hand  side  of  the  unit.  The  trigger 
uses  TTL  logic,  and  is  activated  when  the  signal  level  is  connected  to 

0  ground.  This  made  fabricating  a  trigger  quite  easy.  A  coax  cable  was 

attached  to  the  trigger  input.  A  short  piece  of  insulated  wire,  with 
one  end  stripped  of  insulation,  was  inserted  into  the  female  connector 

0  (signal)  of  the  BNC  fitting  at  the  free  end  of  the  cable.  To  activate 

the  trigger,  the  wire  was  grounded  to  the  shield  of  the  coax  cable. 
(27:5-1,9) . 

• 
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After  the  trigger  event,  the  input  channels  are  recorded  until  the 
buffers  are  full  at  time  t  Each  buffer  recorded  can  then  be 

tot 

plotted.  A  special  function  allows  the  buffer  to  be  stepped  through,  one 
location  at  a  time.  The  plotter  moves  the  pen  directly  over  the  point 
on  the  plot  representing  the  current  storage  location;  its  value  is 
observable  on  the  LCD  display  on  the  plotter.  If  the  point  represents  a 
point  of  interest,  such  as  a  peak  in  a  curve  or  a  zero  amplitude,  a 
function  can  be  selected  to  mark  the  point  with  a  '+'  and  annotate  it 
with  voltage  and  time  values  (27:4-1,10). 


Determination  of  Gain  and  Phase  Angle  Using  the  Model  7090A 
To  interface  input  signals  with  the  model  7090A  use  banana  plugs, 
or  a  coax  cable  with  a  dual  banana  plug  adaptor.  Connect  the  ground  of 
the  device  or  circuit  being  measured  to  the  ground  jack  on  the 
recorder.  Configure  the  input  signal  from  the  function  generator  into 
channel  1  of  the  recorder,  and  the  output  of  the  Oldham-Zoski  circuit 
into  channel  2.  Repeat  the  constant  amplitude  input  voltage  check  as 
previously  mentioned.  Input  values  can  be  observed  on  the  LCD  display 
and  do  not  need  to  be  plotted  (this  was  repeated  once;  the  voltage 
remained  constant  over  the  frequency  range) . 

For  a  given  frequency  determine  the  total  time  from  Eq  (FI)  and  the 
desired  pre-trigger  time  and  input  to  the  plotter.  Select  channels  1 
and  2  for  recording.  Set  the  expected  voltage  for  each  channel.  Set 
the  frequency  generator  to  the  desired  value  and  verify  with  the 
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counter  timer.  Turn  the  frequency  generator  off.  This  permits 
recording  of  the  start-up  transient  of  the  fractional  differentiator. 

Press  the  ‘FILL  BUFFER"  switch;  after  sufficient  time  has  passed  to  load 
the  pre-trigger  buffer,  turn  on  the  function  generator  and  ground  the 
trigger  wire.  After  the  buffers  fill,  plot  the  buffers  twice 
(27:4-7,8).  On  one  annotate  the  points  of  interest  on  the  input 
channel.  On  the  other  annotate  the  points  of  interest  for  the  circuit 
output.  The  points  of  interest  used  for  the  data  in  Appendix  G  are  the 
maximum  and  minimum  points  of  each  cycle  and  the  negative  slope  time  axis 
crossings  (the  LCD  display,  in  the  annotation  mode,  can  be  invaluable  in 
selecting  these  points)  (27:4-8). 

The  period,  gain,  and  phase  are  calculated  as  follows: 

Period: 

FINAL  PEAK  VOLTAGE  TIME  <  SEC  >  -  INITIAL  PEAK  VOLTAGE  TIME  CSJEC) 

NUMBER  OF  CYCLES  BETWEEN  THE  PEAK  VALUES  USED  (F4) 

Gain : 

Z  [output  voltage  peaks  >|qutput  VOLTAOE  MINIMUMS  |]/  [tf  MAXES  ♦  mins] 

Z  [INPUT  VOLTAOE  PEAKS  +  |  INPUT  VOLTAOE  MINIMUMS  |  ]  /  MAXES  +  MINs] 

(F5) 

This  method  compensates  for  offsets  in  the  input  and  output  voltages. 

Phase  Angle  (using  paired  input  and  output  time  axis  crossings) : 

Z  [INPUT  TIME  AXIS  CROSSINGS  (a«C)  7  -  Z  [  OUTPUT  TIME  AXIS  CROSSINGS] 

NUMBER  OF  CROSSINGS  EVALUATED  *  TOTAL  PERIOD  (F6) 
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This  method  proved  much  easier  and  faster  than  observing  the 
oscilloscope,  counter/timer ,  and  continually  making  the  necessary 
adjustments . 

Analog  Components  and  Tasks 

Figure  20  shows  the  PACE  analog  computer  (model  TR48) ,  with  its 
three  major  panel  subdivisions  -  the  control  and  monitoring  panel  on  the 
left,  the  component  panel  in  the  center,  and  the  potentiometer 
adjustments  on  the  right. 

Control  and  Monitoring  Panel  (22:1-11) 

There  are  three  items  of  concern  for  operating  the  analog  computer 
on  this  panel;  the  multi-colored  mode  select  keys,  digital  voltmeter 
(DVM) ,  and  the  operational  ampl i f ier/potentiometer  selector  keypad. 

First  the  mode  select  keys.  There  are  four  of  interest  for  the  tasks  at 
hand;  'PS‘,  potentiometer  set;  *IC',  initial  condition;  "OP',  operation, 
and  "HD",  hold  modes.  To  select  a  specific  mode,  simply  depress  the 
appropriate  key.  The  amplifier/potentiometer  selector  keypad  simply 
selects  one  of  the  numbered  potentiometers  or  amplifiers  (depending  on 
whether  the  'A'  or  'P'  selector  is  chosen)  for  display  on  the  DVM. 

To  view  the  value  of  a  specific  potentiometer  used  for  establishing 
the  coefficients  of  the  differential  equations,  depress  the  ”PS'  key. 
This  isolates  the  potentiometers  from  their  input  voltages,  giving  their 
true  values  in  the  circuit.  Depress  the  black  "P‘  key,  then  the 
appropriate  two-digit  potentiometer  identifier  as  listed  on  the 
component  panel.  The  desired  value  can  now  be  read  off  the  DVM. 


(159) 


To  view  the  output  value  of  a  specific  amplifier  in  any  mode, 
depress  the  ’A"  key  and  appropriate  two  digit  designator  of  the  desired 
amplifier  as  found  on  the  component  designator  panel.  The  output  of  all 
amplifiers  is  a  varying  voltage  level.  The  DVM  represents  this  level 
as  a  certain  percentage  of  the  computer’s  reference  voltage  of  10  volts. 
Therefore,  when  the  DVM  reads  0.2500,  the  output  of  a  given  amplifier  is 
actually  2.5  volts  (22:1-11) 

When  are  the  different  modes  used’  Use  'PS'  mode  when  setting  and 
rechecking  potentiometer  values.  The  potentiometers  have  a  tendency  to 
drift,  so  after  a  few  data  runs  switch  back  into  the  'PS'  mode  and 
recheck  the  values.  Use  *IC”  at  the  completion  of  a  data  run,  to  abort  a 
run,  or  to  reestablish  initial  conditions.  When  in  the  'IC'  mode,  the 
initial  conditions  for  each  operational  amplifier  and  potentiometer  can 
be  viewed  on  the  DVM.  This  capability  is  ideal  for  verifying  proper 
circuit  configuration.  After  everything  is  properly  configured, 
depressing  the  "OP*  button  starts  the  simulation  by  activating 
the  integrators.  The  “HD”  key  freezes  all  computations  at  the  current 
levels  and  permits  observation  of  intermediate  values  (22:1-20). 

Located  on  the  center  panel  of  the  model  TR48  are  many  components 
required  to  simulate  control  systems  and  differential  equations.  Only 
the  ones  used  for  the  open-  and  closed-loop  simulations  will  be 


described . 


These  functions  are  provided  by  the  '6.514  Dual  DC  Amplifier' 
(22:3-1).  All  amplifiers  on  the  model  TR48  are  inverting  amplifiers, 

•  that  is,  the  output  voltage  is  the  negative  of  the  input  voltage.  A 

voltage  input  to  the  jack  labeled  ’10"  produces  a  negative  ten  times  the 
input  voltage.  If  more  than  one  voltage  is  input  to  the  operational 
p  amplifier,  the  output  is  the  sum  of  the  input  voltages  multiplied  by  the 

value  of  the  input  jacks  used. 

Attenuators  (Potentiometers) 

Without  attenuators,  precision  multiplication  would  not  be 
possible.  The  "42.283  Attenuator  Module"  (22:2-1)  consists  of 
5,  wire-wound  resistors  and  whose  resistance  is  established  by  the 
setting  associated  with  its  control  dial.  The  dials  are  divided  into 
one  ten- thousandths  increments.  Circuit  schematics  and  flowchart 
symbols  for  the  potentiometers  are  shown  in  Figure  22.  To  obtain  a 
desired  voltage,  for  example  0.8973  x  the  following  procedure  applies: 

1.  Patch  the  input  voltage  into  the  upper  terminal  of  the 
potentiometer 

2.  Patch  the  output  of  the  potentiometer  to  the  input  of  the  next 
device 

3.  Determine  the  component  identifier  of  the  potentiometer  (for 
example  P07) 


4.  Depress  the  'PS'  key  to  go  into  Potentiometer  Set  mode 

5.  Select  'P',  "0",and  '7',  using  the  selection  keys 


6.  Observe  the  current  value  of  the  pot  on  the  DVM 

7.  Unlock  the  'P07'  dial  on  the  right  panel  and  adjust  until 


the  DVM  reads  0.8973 

8.  Lock  the  dial,  making  sure  the  potentiometer  value  remains  the 
same 

The  output  of  P07  will  now  be  0.8973  of  the  value  input. 

Integrators 

The  ‘12.1322  Dual  Integrator  Network"  (22:4-6)  produces 
the  integral  of  the  input  voltage  with  respect  to  time.  Reference  the 
symbol  under  ‘X(0)  =2.0' ,  depicted  in  Figure  34,  for  the  circuit 
schematic  associated  with  an  integrator.  Note  the  value  at  the  top  of 
the  x  integrator,  this  is  the  initial  condition  applied  to  the 
integrator,  and  is  the  output  of  the  integrator  when  in  the  "IC"  mode. 
This  value  is  established  with  a  reference  voltage  input  to  the  ‘IC  jack 
on  the  integrator.  Again,  the  integrator  is  an  inverter  and 
the  negative  value  of  this  voltage  is  the  initial  condition.  When  placed 
in  the  ‘0P‘  mode,  the  output  of  the  integrator  is: 

-  (  [  input  voltage(t)  ]  dt  +  IC  (F7) 

Comparator 

To  successfully  interface  the  model  7090A  with  the  analog  computer, 
the  recorder  must  receive  its  input  signal  at  the  moment  the  op  button 
is  depressed  on  the  control  panel.  As  previously  discussed,  the 
model  7090A  activates  its  trigger  when  the  signal  voltage  is  connected 
to  ground.  To  correctly  patch  this  logic,  the  "40.404  Relay  Comparator" 


(shown  schematically  in  the  upper  right  hand  corner  of  Figure  23 
isolating  the  one-half  derivative  input)  is  used.  The  trigger  input 
signal  from  the  model  7090A  recorder  is  patched  into  the  common  terminal 
of  the  switching  relay,  the  trigger  ground  into  the  comparator  terminal 
marked  *-*.  The  switching  voltage  required  to  initiate  the  trigger 
signal  is  provided  at  the  output  of  the  'OP'  bus  jack  on  the 
12.132  Dual  Integrator  when  the  'OP'  mode  is  selected.  When  the  computer 
is  in  ‘IC*  mode  the  output  of  the  ‘OP*  bus  jack  =  0.  This  makes  the 
output  of  the  comparator  positive  and  the  model  7090A  trigger  signal  is 
isolated  from  ground.  When  the  ‘0P‘  mode  is  selected,  the  output  of  the 
comparator  is  negative,  closing  the  connection  between  signal  and 
ground . 

Function  Switches 

On  the  right-hand  panel  of  the  model  TR48  analog  computer,  below  the 
potentiometer  adjustment  dials,  are  five  toggle  switches.  These 
switches  are  three-position,  single-poe  switches  connected  to  the 
‘12.766  Function  Switch  Module’  (22:2-6).  When  the  switch  is  moved  to 
the  left,  the  input  is  connected  to  the  ‘L‘  output  jack.  Likewise,  the 
‘R‘  contact  is  closed  when  the  switch  is  moved  to  the  right.  This 
switch  is  useful  in  connecting  different  forcing  functions  to  a  system 
model.  In  this  specific  case,  several  f(t)s  could  be  available  for 
switching  in  and  out  of  the  circuit.  This  switch  was  used  in  the 
‘Total-Cycle  Simulation"  to  apply  and  remove  the  input  step  voltage. 


More  detailed  information  on  these  components  can  be  found  in  the 


EAI  PACE  (model  TR48)  Analog  Computer  Reference  Handbook  (22) . 

Step3  for  Programming  an  Example  Simulation  (23 : 84+) 

For  an  example  simulation  set-up,  program  a  harmonic  oscillator 
+  x  =  0  with  an  initial  condition  x(0)=2.5  .  To  run  the  simulation: 

1)  Select  three  operational  amplifiers  and  3  potentiometers  and 
define  as  follows: 


Component 
AO  1 


A02 

A03 

P01 

P02 

P03 


\  confi 
J  integ 


gure  as 
tegrators 


Purpose 

x 

x 

X 

m 

k 

x  (0) 


2)  Make  the  following  patches: 


From 

To 

AO  1  Output 

P01 

input 

P01  Output 

xlO 

A02  input 

A02  Output 

xl 

A03  input 

A03  Output 

P02 

input 

10  Ref  on  P00-04  module 

?03 

input 

P03  Output 

A03 

IC 

P02  Output 

x  10 

A01  Input 

A03  Output  HP7G90A  Channel  1  input 

3)  Patch  the  ground  jacks  for  A01,  A02,  and  A03  and  the  model 

7090A  Channel  1  grounds  together  to  eliminate  floating  grounds. 

4)  Configure  the  '40.404  Comparator'  and  model  7090A  per  the 
description  under  comparators  in  this  appendix. 

5)  Depress  ’PS"  selector  on  the  control  console 
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6)  Turn  the  model  7090A  recorder  and  the  model  TR48  analog 


computer  on 

7)  Press  "Restore  Setup'  on  the  model  7090A  to  establish  standard 

values.  Make  desired  changes  to 

Channel  1  Range 
Channel  1  Offset 
Trigger  Value 
Total  Time 

8)  For  each  potentiometer,  press  the  appropriate  component  selection 
keys.  Unlock  the  associated  potentiometer  set  dials,  and  establish  the 
following  values  on  the  DVM  for  each  potentiometer: 


Potentiometer 

Value 

P01 

0.  1000 

P02 

0. 1000 

P03 

0.2500 

9)  Depress  the  "IC"  mode  select  button. 

10)  Select  "A03‘  for  display  on  the  DVM.  It  should  read  =  0.2500 
(whatever  the  setting  on  P03  was) 

11)  Press  "BUFFER  FILL"  on  the  model  7090A,  and  allow  adequate 
time  for  the  pre-trigger  buffer  to  fill 

12)  Depress  "OP"  on  the  model  TR48  analog  computer 

13)  When  the  "Buffer  Full"  light  comes  on  on  the  model  7090A: 

a)  Press  "IC"  on  the  model  TR48  to  end  the  run 

b)  Insert  paper  in  the  model  7090A 


> 
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c)  If  a  grid  is  desired  on  the  plot 


1)  Select  a  pen  color  with  the  ‘PEN  SELECT'  button 

2)  Press  'Grid' 

®  14)  When  plotting  is  complete,  annotate  the  points  of  interest  as 

previously  mentioned. 

The  resulting  output  should  be  a  cosine  function  with  amplitude  2.5 

%  volts  and  period  (jt/2)  Hz.  This  circuit  can  be  used  for  the 

analog/half-order  compatibility  verification  trial.  Similar  methods  can 
be  used  for  simulating  other  second-order  equations.  In  this  simulation 

•  P01  and  P02  were  not  required  to  obtain  the  coefficients  m  and  k  since 
they  were  both  identically  equal  to  1.  The  output  of  A01  and  A03  could 
have  been  connected  directly  to  the  next  component  with  the  same 

•  results.  However,  the  example  provided  an  exercise  for  establishing  the 
proper  coefficients  for  each  term. 
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Oldaham  -  Zoski  Circui t  Val idation  Data 


Appendix  G 

Oldaham  2.  Zoski  Circui t  Val idation  Data 
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I 

I 


Table  4.  Circuit  1  Performance  Validation  Data 


Frequency  (Hz) 

Gain  (dB) 

Phase  (deg) 

0.0116 

-47.85 

44.74 

0.0505 

-41.71 

45.81 

0.0963 

-38.95 

45.76 

0.1451 

-37.26 

47.00 

0.1934 

-36.03 

48.04 

0.5102 

-31.67 

46.53 

0.9780 

-28.88 

45.77 

5.042 

-22.08 

44.16 

9.750 

-19.36 

45.40 

Table  5.  Circuit  2  Performance  Validation  Data 


Frequency  (Hz) 

Gain  (dB) 

Phase  (deg) 

0.0107 

-38.96 

42.18 

0.0503 

-32.75 

50.35 

0.0976 

-29.34 

48.51 

0.1470 

-27.61 

46.97 

0. 1944 

-26.62 

44.08 

0.5034 

-23.53 

40.11 

0.9900 

-21.03 

43.20 

5.012 

-13.89 

46.92 

9.670 

-10.79 

46.  19 
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Figure  48. 
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Figure  49.  Cirt 
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CIRCUIT  1  -  FREQUENCY  =  .0963  HZ 


m  VO  ITS 
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-460 


CIRCUIT 


CIRCUIT  1  -  FREQUENCY  =  .1451  HZ 


Figure 55 .  Circuit  1  Validation  TIME  (SEC) 


CIRCUIT  1  -  FREQUENCY  =  .1934  HZ 


CIRCUIT 


Figure  58.  Circuit  I  Validation  ^  (SEC) 
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CIRCUIT  1  -  FREQUENCY  =  9.750  HZ 


(185) 


Validation  "^E  (SEC) 


CIRCUIT 


Figure  67.  Circuit  2  Validation  ^ 1 M  ^  ^  ^  ^ ) 


Figure  66.  Circui 1 2  Validation  TIME  (SEC) 
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Figure  77.  Czra/i  £  Validation  ^  ^  E  C ) 
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Figure  78 .  Circuit  2  Validation  (SEC) 


CIRCUIT  2  -  FREQUENCY  =  9.670  HZ 


Figure  82.  Circuit  2  Validation 


CIRCUIT  2  -  FREQUENCY  =  9.670  HZ 


Figure  83.  Circuit  2  Validation 
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an  input  signal  are  investigated.  The  method  selected  employs  a  circuit 
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integration.  The  results  presented  show  the  circuit  accurately 
differentiates  a  sinusoidal  input  for  a  frequency  range  spanning  0.01  Hz 
to  10.0  Hz. 

The  second-order  differential  equation  above  is  simulated  on  an 
analog  computer.  An  optimal  u(t)  is  then  used  for  feedback  modification 
of  the  original  open-loop  system.  Improved  system  performance  resulted. 

A  Laplace  transform  method  and  a  Mittag-Lef f ler  expansion  provide 
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